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ABSTRACT 

We present an implementation of stellar evolution and chemical feedback for smoothed parti- 
cle hydrodynamics (SPH) simulations. We consider the timed release of individual elements 
by both massive (Type II supernovae and stellar winds) and intermediate mass stars (Type la 
^ c"| supernovae and asymptotic giant branch stars). We illustrate the results of our method using 

a suite of cosmological simulations that include new prescriptions for radiative cooling, star 
formation, and galactic winds. Radiative cooling is implemented element-by-element, in the 
presence of an ionizing radiation background, and we track all 1 1 elements that contribute 
significantly to the radiative cooling. 

While all simulations presented here use a single set of physical parameters, we take 
specific care to investigate the robustness of the predictions of chemodynamical simulations 
with respect to the ingredients, the methods, and the numerical convergence. A comparison of 
nucleosynthetic yields taken from the literature indicates that relative abundance ratios may 
only be reliable at the factor of two level, even for a fixed initial mass function. Abundances 
relative to iron are even more uncertain because the rate of supernovae la is not well known. 
[ We contrast two reasonable definitions of the metallicity of a resolution element and find 

that while they agree for high metallicities, there are large differences at low metallicities. 
We argue the discrepancy is indicative of the lack of metal mixing caused by the fact that 
metals are stuck to particles. We argue that since this is a (numerical) sampling problem, 
solving it using a poorly constrained physical process such as diffusion could have undesired 
| consequences. We demonstrate that the two metallicity definitions result in redshift z = 

stellar masses that can differ by up to a factor of two, because of the sensitivity of the cooling 
rates to the elemental abundances. 

Finally, we use several 512 3 particle simulations to investigate the evolution of the distri- 
• bution of heavy elements, which we find to be in reasonably good agreement with available 

observational constraints. We find that by z = most of the metals are locked up in stars. The 
gaseous metals are distributed over a very wide range of gas densities and temperatures. The 
shock-heated warm-hot intergalactic medium has a relatively high metallicity of ~ 10 _1 Zq 
that evolves only weakly and is therefore an important reservoir of metals. Any census aiming 
to account for most of the metal mass will have to take a wide variety of objects and structures 
into account. 

Key words: cosmology: theory — galaxies: abundances — galaxies: formation — intergalac- 
tic medium — methods: numerical 
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1 INTRODUCTION 

Nucleosynthetic processes within stars and supernovae (SNe) 
change the abundances of elements in the Universe over time. 
* E-mail: wiersma@strw.leidenuniv.nl As stars release these elements back into the diffuse interstellar 
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medium (ISM) that surrounds them, subsequent generations of stars 
are born with different chemical compositions. On larger scales, 
winds driven out of galaxies by the energy released by dying stars 
and/or active galactic nuclei modify the composition of the inter- 
galactic medium (IGM) and hence that of future galaxies. The IGM 
within groups and clusters of galaxies exerts a ram-pressure on their 
member galaxies as they move along their orbits. This head-wind 
can be sufficiently strong to strip the ISM from the galaxies, thus 
mixing the elements that were present in the ISM into the IGM. 
These and other processes make the (re-)cycling of elements over 
time through stars, galaxies and diffuse gas a complicated problem 
that is well-suited for studies using hydrodynamical simulations. 

Simulating the exchange of elements between stars, galaxies, 
and their environments is the subject of this paper. Specifically, this 
paper describes and tests a method to follow the timed release and 
subsequent spreading of elements by stellar populations formed in 
cosmological, hydrodynamical simulations. 

The elemental abundances of stars, galaxies, and diffuse gas 
are of great interest for a variety of reasons. First, the dissipation 
of binding energy by the emission of "cooling" radiation is a fun- 
damental ingredient of models of the formation of both stars and 
galaxies. Because radiative cooling rates are very sensitive to abun- 
dance variations, the rate at which gas can collapse into galaxies 
is to a large extent determined by its chemical composition. Simi- 
larly, the initial mass function (IMF) of stars may well depend on 
the chemical composition of the gas from which they formed. 

Second, knowledge of elemental abundances can give great in- 
sight into a variety of key astrophysical processes. For example, the 
distribution of elements in the IGM provides stringent constraints 
on models of galactic winds and ram-pressure stripping in clusters. 
The relative abundances of elements can tell us about the IMF of 
the stars that produced them. The abundance of alpha elements - 
which are released on short time-scales through core collapse SNe 
- relative to iron - most of which is released much later when in- 
termediate mass stars explode as type la SNe - tells us about the 
time-scale on which stars have been formed. On the other hand, the 
absolute abundances of metals provide us with a measure of the 
integral of past star formation. 

Third, emission and absorption lines of individual elements 
constitute key diagnostics of physical conditions, such as gas den- 
sities, temperatures, radiation fields and dust content. The strengths 
of said lines depend on a number of factors, but invariably the abun- 
dances play an important role. 

Fourth, individual lines are the only observable signatures for 
most of the diffuse gas in the Universe. For example, oxygen lines 
are the most promising tools to detect the warm-hot IGM that may 
host a large fraction of the cosmic baryons. To verify this prediction 
of structure formation in a cold dark matter universe, it is necessary 
to know the abundance of oxygen in this elusive gas phase. 

Clearly, following the timed release of the elements by stars 
and their subsequent dispersal through space, a collection of pro- 
cesses sometimes referred to as "chemodynamics", is a critical in- 
gredient of any realistic simulation of the formation and evolution 
of galaxies. The implementation of chemodynamics into cosmo- 
logical simulations is a challenging problem, mainly because such 
models lack the resolution to resolve many of the relevant physical 
processes and hence require "sub-grid" recipes. 

The first study to implement metal enrichment into a 
Smoothed Particle Hydrodynamics (SPH) code did not distin- 
guish between different elements and only considered enrichment 
by c ore collapse SNe in the ins tantaneous recycling approxima- 
tion dSteinmetz & Muellej[l994l) . Since then many authors have 



implemented mo re sophisticated recipes for chemodynamics into 
SPH c odes (e. 

1200 ll: iKawatal 120011 ; iLia et alj 120021 : iKawata & Gibson! 12003 



.. .20011: Iliaet alJ 120021: iKawat a & Gibson 
Valdar ninil 120031; iKobavashil |2004 ISommer-Larsen et alj |2005 



iTornatore et al J 120071 : lOppenheimer & Pavel I2008I) . ranging from 
recipes that distinguish only between SN type la and SN type 
II elements to models that follow large numbers of individ- 
ual elements released by AGB stars, SNe la, SNe II, and the 
winds from their progenitors. Three-dimensional chemodynami- 



of individual objects such as galaxies (e.f 


:. TheisetalJJl992; 


Raiteri et al.ll9961;lBerczikll999: Recchi et al 


120011: Kawatal2001 ; 


Kawata & Gibson 2003l:lKobavashill2004; Mori & Umemura 2006; 


Govemato et al. 20070 and clusters of galaxies (e.g. |Lia et al. 
2002; Valdarnini 2003; Tornatore et al. 2004; Sommer-Larsen et al. 


2005 
2008 


; Romeo et al. 2006; Tornatore et alj|2007l; see Borgani et al. 


for a review), but 


in recent years chemodynamical cos- 



mological simulation s have also become mo r e common (e.g. 



i Mosconi et alj 120 01 ; Scannapi eco et alj 1 20051 ; iKobavashi et alj 
l2007i: lOppenh eimer & Dave 2008). Note t hat several recent cos- 
mological works 1 Scan napieco et al. 2005; Kobavashi et alj|2007l : 
Oppenheim er & Pavel 2008 ) use our hydrodynamica l code of 
choice, namely the SPH code GADGET jspringel 2005), although 
they use different implementations of chemodynamics, star forma- 
tion, galactic winds, and radiative cooling. 

Here we present a new implementation of chemodynamics 
that follows the timed release - by AGB stars, SNe la, SNe II, 
and the winds from their progenitors - of the 1 1 elements that 
IWiersmaet all (2009) found to contribute significantly to the radia- 
tive cooling at T > 10 4 K. We present a subset of high-resolution 
cosmological, hydrodynamical simulations taken from the Over- 
whelmingly Large Simulations (OWLS) project (Schaye et al., in 
prep.) and use them to test the robustness of our prescription to 
numerical parameters and to make some predictions regarding the 
large-scale distribution of heavy elements. These simulations in- 
clude element-by-element cooling, a first for cosmological chemo- 
dynamical simulations, and they take the effect of photo-ionization 
by the UV/X-ray background radiation on the cooling rates into 
account, which is normally either ignored or included only for hy- 
drogen and helium. While this paper focuses on a single implemen- 
tation of the relevant physics, we will present the effects of varying 
the physical parameters (e.g., cosmology, star formation, cooling, 
feedback from star formation) in a future paper that will make use 
of a larger fraction of the OWLS data set. 

This paper is organized as follows. The simulations and the 
prescriptions used for star formation, feedback from star forma- 
tion, and radiative cooling are summarized in Sj2] In Sj3]we discuss 
in some detail the ingredients that we take from models of stel- 
lar evolution or observations, including the IMF, stellar lifetimes, 
stellar yields, and SN la rates. In Sj4]we discuss how we combine 
all these ingredients to predict the mass released by a simple stellar 
population as a function of its age and present some results. Our im- 
plementation of chemodynamics into SPH is discussed in Sj5] This 
section also contrasts two possible definitions of elemental abun- 
dances in SPH and demonstrates that the choice of definition can 
be extremely important due to limitations intrinsic to SPH. In ij6] 
we use our simulations to address the question "Where are the met- 
als?" and we show how the answer is influenced by the definition 
of metallicity that is used. Finally, we provide a detailed summary 
in Sj7]and investigate convergence with the size of the simulation 
box and with resolution in appendices Q5] and [C] respectively. 

For the solar abundance we use the metal mass fraction 
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Table 1. Adopted solar abundances. 



Element 






Mass Fraction 


H 


1 




0.7065 


He 


0.1 




0.2806 


C 


2.46 x 10" 


-4 


2.07 x 10 -3 


N 


8.51 x 10" 


_5 


8.36 x 10~ 4 


O 


4.90 x 10" 


-4 


5.49 x 10~ 3 


Ne 


1.00 x 10" 


-4 


1.41 x 10~ 3 


Mg 


3.47 x 10" 


-5 


5.91 x 10~ 4 


Si 


3.47 x 10" 


-5 


6.83 x 10" 4 


S 


1.86 x 10" 


-5 


4.09 x 10" 4 


Ca 


2.29 x 10" 


-6 


6.44 x 10~ 5 


Fe 


2.82 x 10" 


-5 


1.10 x 10~ 3 



Zq = 0.0127, corresponding to the value obtained using the 
de fault abundance set of CLOUDY (version 07.02; last described 
by iFerland et alj Il998l) . This abundance set ( see T able Q} com- 
bines t he abundances of lAllende Prieto et alj hoOll |2002) for C 
and O, Holwegerj d200ll) for N, Ne, Mg, Si, and Fe, and assumes 
nHe/n-H = 0.1 which is a typical value for nebular with near-solar 
compositions. Note that much of the literature assumes Zq — 0.02, 
which is 0.2 dex higher than the value used here. 



2 SIMULATIONS 

We performed cosmological, gas-dynamical simulations using a 
modified version of the iV-body Tree-PM, SPH c ode GADGET II I, 
which is a modified version of the code GADGET II ( Springel 2005). 
Our prescriptions for star formation, for feedback from core col- 
lapse supernovae and for radiative cooling and heating are sum- 
marized below. Our prescription for stellar evolution, which is the 
subject of this paper, is described in detail in the next section. 

We use the suite of simulations of varying box sizes and parti- 
cle numbers listed in Table|2]together with the corresponding parti- 
cle masses and gravitational force softening scales. We use fixed 
comoving softening scales down to 2 = 2.91 below which we 
switch to a fixed proper scale. Our largest simulations use 2 x 512 3 
particles in boxes of comoving size L = 25, 50, and 100 ft" 1 Mpc. 

The i nitial particle positions and velocities are obtained from 
glass-like I White 1 19941) initial condit ions using CMBFAST (ver- 
sion 4.1; ISeliak & Zaldarri aga 1996) and employing the Zel- 
dovich approximation to linearly evolve the particles down to 
redshift z = 127. We assume a flat ACDM universe and em- 
ploy the set of cosmological parameters [£l m ,os,n 3 ,h] — 
[0.238,0.0418,0.762,0.74,0.951,0.73], as determined from the 
WMAP 3-year data a nd consistent with the WMAP 5-year data 
(Komatsuetal. 2008). In addition, the primordial helium mass 
fraction was set to Ah = 0.248. 

We employ t he star formation recipe of 
ISchave & Dalla Vecchij (2008), to which we refer the reader 
for details. Briefly, gas with densities exceeding the critical 
density for the onset of the thermo-gravitational instability 
(n H ~ 10~ 2 - IP" 1 cm ) is expected to be multiphase and 
star-forming dSchavd [2004h . We therefore impose an effective 
equation of state with pressure P oc p% eH for densities exceeding 
71h =0.1 cm -3 , normalized to P/k = 10 3 cm" 3 K for an atomic 
gas at the threshold density. We use y e g = 4/3 for which both the 



Our value of erg is 1.6 cr lower than allowed by the WMAP 5-year data. 



Jeans mass and the ratio of the Jeans length and the SPH kernel are 
independent of the density, thus preventing spurious fragmentation 
due to a lack of numerical resolution. 

The local Kennicutt-Schmidt star formation law is analytically 
c onverted and implemented as a p ressure law. As we demonstrated 
in lSchave & Dalla Vecchial (2008). our method allows us to repro- 
duce arbitrary input star formation laws for any equati on of state 
withou t tuning any parameters. We use the observed Kennicutt 
(1998) law, 



S, = 1.5 x lO" 4 M0yr _1 kpc" 



1 M pc" 2 



(1) 



where we divided Kennicutt's normalization by a factor 1.65 to ac- 
count f or the fact that it assumes a Salpeter IMF whereas we are 
using a lChabried J2003I) IMF. 

Galactic winds a re im plemented as described in 
iDalla Vecchia & Schavd (2008). Briefly, after a short delay 
of isN = 3 x 10 7 yr, corresponding to the maximum lifetime 
of stars that end their lives as core-collapse supernovae, newly- 
formed star particles inject kinetic energy into their surroundings 
by kicking a fraction of their neighbours in a random direction. 
The simulations presen te d her e use the default parameters of 
IDalla Vecchia & Schavel d2008h . which means that each SPH 
neighbour i of a newly-formed star particle j has a probability 
of Tfrrij/ y^ J ^=s b mi of receiving a kick with a velocity v w . We 
choose r\ = 2 and ?j w = 600 kms (i.e., if all baryonic particles 
had equal mass, each newly formed star particle would kick, on 
average, two of its neighbours). Assuming that each star with 
initial mass in the range 6 — 100 Mq injects 10 51 erg of kinetic 
energy, these parameters imply that the total wind energy accounts 
for 40 per cent of the available kinetic energy for a Chabrier 
IMF and a stellar mass range 0.1 — 100 Mq (if we consider 
only stars in the mass range 8 — 100 Mq for type II SNe, this 
works out to be 60 per cent). The value 77 = 2 was chosen to 
roughly reproduce the peak in the cosmic star formation rate. 
Note that contrary to the w idely-used kinetic feedback recipe of 
ISpringel & Hernquistl J2003l) . the kinetic energy is injected locally 
and the wind particles are not decoup l ed h ydrodynamically. As 
discussed by IDalla Vecchia & Schav these differences 

have important consequences. 

Radiative cooli ng was implemented according to 
IWiersma et ail [2009f In brief, net radiative cooling rates 
are computed element-by-element in th e presence of the c osmic 
microwave background (CMB) and a Illaardt & M adaul d200l[ 
hereafter HM01) model for the UV/X-ray background radiation 
from quasars and galaxies. The contributions of the eleven elements 
hydrogen, helium, carbon, nitrogen, oxygen, neon, magnesium, 
silicon, sulphur, calcium, and iron are interpolated as a function of 
density, temperature, and redshift from tables that have been pre- 
computed using the publi cly available photo-i onization package 
CLOUDY, last described bv lFerland et al. ( 1998), assuming the gas 
to be optically thin and in (photo-)ionization equilibrium. 

Hydrogen reionization is implemented by turning on the 
evolving, uniform ionizing background at redshift z = 9. Prior to 
this redshift the cooling rates are computed using the CMB and a 
photo-dissociating background which we obtain by cutting off the 
z = 9 HM01 spectrum at 1 Ryd. Note that the presence of a photo- 
dissociating background suppresses H2 cooling at all redshifts. 



We used their equation (3) rather than (4) and CLOUDY version 05.07 
rather than 07.02. 
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Table 2. List of simulations used. The columns give the comoving size of the box L, the total number of particles per component N (dark matter and 
baryons), the initial baryonic particle mass mj,, the dark matter particle mass r«dm. the (Plummer-equivalent) comoving softening length e C om, the maximum 
(Plummer-equivalent) proper softening length e prop , and the redshift at which the simulation was stopped z end . 



Simulation 


L 


TV 




"Mm 


e com 


e prop 


^cnd 




(h- 1 Mpc) 




(Hm ) 


(/i- 1 M Q ) 


(ft- 1 kpc) 


(h- 1 kpc) 




L006N128 


6.25 


128 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.5 


2 


L012N256 


12.50 


256 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.5 


2 


L025N128 


25.00 


128 3 


8.7 x 10 7 


4.1 x 10 s 


7.81 


2.0 





L025N256 


25.00 


256 3 


1.1 x 10 7 


5.1 x 10 7 


3.91 


1.00 


2 


L025N512 


25.00 


512 3 


1.4 x 10 6 


6.3 x 10 6 


1.95 


0.5 


1 


L050N256 


50.00 


256 3 


8.7 x 10 7 


4.1 x 10 s 


7.81 


2.0 





L050N512 


50.00 


512 3 


1.1 x 10 7 


5.1 x 10 7 


3.91 


1.0 





L100N128 


100.00 


128 3 


5.5 x 10 9 


2.6 x 10 10 


31.25 


8.0 





L100N256 


100.00 


256 3 


6.9 x 10 8 


3.2 x 10 9 


15.62 


4.0 





L100N512 


100.00 


512 3 


8.7 x 10 7 


4.1 x 10 s 


7.81 


2.0 






1 Q 5 ^ , , , , , , , , , , , , 

1 He reionization heat 

No extra heat 




Redshift 

Figure 1. A comparison between the simulated and observed evolution of 
the temperature of gas with density equal to the cosmic mean. The dashed 
curve indicates the temperature evolution for a gas parcel undergoing Hub- 
ble expansion, computed using the same radiative cooling and heating rates 
as are used in our simulations for gas of primordial composition. These net 
cooling rates include the effects of the CMB and, below z = 9, of the 
evolving UV/X-ray background (UVB) from galaxies and quasars as mod- 
eled by HM01, assuming the gas to be optically thin and in ionization equi- 
librium. The redshift of hydrogen reionization was set to z = 9 by turning 
on the UVB at this redshift, resulting in a steep increase of th e temperature 
to T ~ 10 4 K. A comparison with with the data points of Scha ve et alj 
(2000), which were derived from the widths of Lyce absorption lines, shows 
that the simulation underestimates the temperature at z ~ 3, the redshift 
around which the reionization of helium is thought to have ended. To mimic 
the expected extra heat input, relative to the optically thin limit, during he- 
lium reionization, the solid curve shows the result of injecting an extra 2 eV 
per proton, smoothed with a a(z) = 0.5 Gaussian centered at z = 3.5. The 
reduced x 2 f° r the solid and dashed curves is 1.1 and 4.5, respectively. 



Our assumption that the gas is optically thin may lead to an un- 
dere stimate of the temperatu re of the IGM shortly after reionization 
(e.g. lAbel & Haehneltjl999h . Moreover, it is well known that with- 
out extra heat input, hydrodynamical simulations underestimate the 
temperature of the IGM at z > 3, the redshi ft around which he- 
lium r eionization is thought to have ended (e. g .|T heuns et al. 1998, 
1 19991 : iBrvan et alj 1 19991 : 1 Schave et al.l l2000l : iRicotti et al.l l2000h . 
We therefore inject 2 eV per proton, smoothed with a a(z) = 0.5 
Gaussian centered at z — 3.5. Figure [T] shows that with this extra 



energy, the predicted tem perature at the me an density agrees well 
with the measurements of lSchave et alj J200C|> . 



3 INGREDIENTS FROM STELLAR EVOLUTION 

Cosmological simulations cannot at present resolve individual 
stars. Instead, a single particle is taken to represent a population 
of stars of a single age, a so-called simple stellar population (SSP), 
with some assumed stellar initial mass function (IMF). The task 
of stellar evolution modules in cosmological hydrodynamical sim- 
ulations is to implement the timed release of kinetic energy and 
mass by star particles. The feedback processes originating from 
stellar evolution that we will consider are winds from asymptotic 
giant branch (AGB) stars, type la supernovae (SNe), type II (i.e. 
core collapse) SN^fl and the winds from their progenitors. In this 
section, we briefly present the ingredients that go into our model, 
while deferring a discussion of alternatives to appendix lAl 

Although many formulations of the stellar initial mass func- 
tion (IMF) exist in the literature, there is consensus that a simple 
power law pr edicts too many l ow mass stars below 1 Mq. We have 
decided on a IChabrien 120031) IMF - because it is an example of 
an IMF with a low mass turnover that fits the observations within 
the uncertainties. For the parameters we assume M c = 0.079, 
a = 0.69 (see appendix I Alb with mass limits of 0.1 Mq and 
100 M©. While this is a sensible choice, uncertainties in the IMF 
can have a significant effect on the global metal production; 

F or s tellar yields, we cho ose the complete set of iMarigol 
d200lh and lPortinari et alj dl998h. along with th e type la SN yields 
of the W7 model of Thieleman n et alj d2003h . This covers mass 
loss fr om intermedia te mass stars as they pass through the AGB 
phase (Marigj |200jh . which with this yield set has an upper limit 
of 5 M0, as well as mass loss and ex plosive nucleosynthes is that 
occur in high mass stars up t o 100 Mp JPortinari" etal.1 19 9811. The 
yields of Thielemann et al. ( 2003) cover the regime of intermedi- 
ate mass stars in close binaries which end their life as type la SNe. 
These two modes of SNe inject energy into the ISM. We inject en- 
ergy from core collapse SNe in kinetic form as discussed in section 

3 Type lb and Ic SN are also core collapse supernova, but we follow other 
authors and use type II to refer to the entire class of core collapse SN. 

4 The high mass stellar yields were modified following the advice of L. 
Portinari - see appendix I A3 I for more details. 
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[2] and in more detail in lDalla Vecchia &Schavd (!2008). The energy 
from SN la is injected thermally among the SPH neighbors of the 
star particles. 

For a more detailed discussion of the yields we refer the reader 
to appendix [A] where we show, by comparing different sets of nu- 
cleosynthetic yields taken from the literature, that they represent 
an important source of uncertainty. In particular, for AGB stars, as 
well as for the elements heavier than nitrogen in massive stars, the 
yields vary by factors of a few. We make the choice that most suits 
our needs, while keeping in mind the resulting uncertainty. 

As in most recent chemodynamical implementations, we relax 
the instantaneous recycling approximation. This means we apply a 
finite lifetime for each star as a function of mass. In consistency 
with our yield sets, w e use the metallicity-dependant lifetimes of 
IPortinari et alj ( fl998h . As shown in appendix I A2l however. there is 
little difference between published lifetimes. 

Since type la SNe are thought to be a product of binary evo- 
lution, there is no simple 'lifetime' to assign to their progenitors. 
Instead, theoretically or empirically derived rates must be used. 
We apply a modification of a e-folding delay function that is nor- 
malised to the cosmic type la SNe rate. In appendix I A4I we show 
that this is yet another source of uncertainty in chemical evolution 
models. Since the precise progenitors are still unknown, and since 
the observations of type la SNe are inconsistent, there is freedom 
to change the total number of type la SNe by factors of a few. 



4 THE MASS EJECTED BY A SIMPLE STELLAR 
POPULATION 

The mass that is ejected during a single time step (t, t + At) by a 
star particle of initial mass m*,o and metallicity Z that was created 
at time t* < t is given by 

Am* = m»(t) -m»(t + At) (2) 

nM z (t-t,) 

= m», / $>(M)m e i(M,Z)dM, (3) 

JM z (i-l,+A() 

where Mz(j) is the inverse of the lifetime function tz(M), 
m e j(M, Z) is the mass ejected by a single star and $(M) is the 
IMF, normalized such that J M$(M)dM = 1. It is straight- 
forward to generalize the above equations to give Am.j, the 
mass ejected in the form of element j. The functions Tz(M) and 
mj, e j (M, Z) need to be taken from stellar evolution and nucle- 
osynthesis calculations provided in the literature. 

4.1 Implementation 

Since stellar evolution and nucleosynthesis depend on initial com- 
position, a correct treatment of chemodynamics would require 
yields (i.e. the ejected mass of various elements) to be tabulated in 
a space with dimensionality equal to the number of elements. Most 
yields are, however, given only as a function of metallicity, where 
metallicity is the total mass fraction in metals. The models used to 
compute the yields assume implicitly that the relative abundances 
of the heavy elements are initially solar. 

One way to implement the release of elements is to interpolate 
the ejected masses, which are given as a function of metallicity, to 
the metallicity of the star particle and to ignore possible differences 
between the relative abundances assumed by the yield calculation 
and the actual relative abundances of the star particle. However, 
this can have some very undesirable consequences. For example, 



because iron does not participate in the nucleosynthesis of interme- 
diate mass stars, the iron that these stars eject was already present 
when the star was born. If an intermediate mass star has a high 
iron abundance for its metallicity, then simple interpolation of the 
yield tables with respect to metallicity would lead to a large un- 
derestimate of the ejected iron mass. That is, it would result in the 
spurious destruction of iron. 

A more accurate way to implement chemodynamics would 
thus be to distinguish metals that were produced/destroyed from 
those that pass through the star without participating in the nucle- 
osynthesis. If we knew, for each ejected element, what fraction of 
the ejected mass resulted from net nucleosynthetic production (i.e., 
production minus destruction) and what fraction simply passed 
through the star (the two fractions should add up to one), then it 
would make sense to assume that the mass that passed through the 
star is proportional to the star's initial abundance. In that case we 
would still neglect the effect of relative abundances on the nucle- 
osynthesis itself (this can only be done by repeating the nucleosyn- 
thesis calculatio ns or using the som ewhat cumbersome Q-matrix 
formalism, e.g. IPortinari et al.|[l998l) . but we would take into ac- 
count the consequences of varying relative abundances on metals 
that do not participate in the nucleosynthesis. 

We have adopted this strategy - as is often done (e.g . 
Recch i et al. 2001; Liaetal. 2002; Sommer-L arsen et alj 120051 ; 
Romeo et al. [20061 ; iTornatore et alj|2007l) - and have implemented 
it as followsQ For a star of given initial mass and metallicity, we 
can write the mass of element j that is ejected into the surrounding 
medium, the sum of two parts: 

m j,cj = "%',0:ej + rnj,p:cy (4) 

Here mj,o:ej is the mass in element j that would have been ejected 
if no nucleosynthesis had taken place and if the star were well- 
mixed, 

(ITlj 0:star \ / n , c \ 

— (m to t,ej) tbl , (5) 

Wltot,0:star / ■ 
/ sim 

where (mj,o :s tar/mtot,o:star) sim is the initial mass fraction of ele- 
ment j in the star as tracked by the simulation and (m, to t,ej) tb i is 
the total ejected mass according to the yield table. Note that we use 
the subscript "sim" to indicate values predicted by the simulation 
and the subscript "tbl" to indicate values that are published by au- 
thors presenting a particular set of stellar yields. The second term 
of equation I0 represents the ejected mass in element j that was 
produced minus destroyed. Assuming that the star is well-mixed, 
we can write it as 

mj, P :ej = (mj, c j) thl - (m t ot,ej) tbl , (6) 

\_ "Hot,0:star J , , 

where (m :)i o:star/m-tot,0:star) tbl is the initial mass fraction of ele- 
ment j that was used for the nucleosynthesis calculation. This term, 
which is often referred to as the yield of element j, will be negative 
if an element effectively gets destroyed (this is for example the case 
for hydrogen). 

Combining equations (4J, (|5), and we obtain: 



5 Note that this strategy applies directly for our chosen yields since they 
are present ed as what is produced min us what was initially in the star. Some 
yields (e.g. Woosley & Weaver 1995) would have to be modified to imple- 
ment this approach. 
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= ( m j,ej)tbl + 



"Hot,0:Bta 



( m 3,0: 3 tar 
\™Hot,0:sta 



(TJtot,ej)tbl- 



(7) 



Note that this expression behaves as expected in the limit that the 
simulated abundances agree with the ones assumed by the yield 
tables (i.e. the second term vanishes). Note also that although the 
abundances appearing in equation ^ are absolute, it is the abun- 
dance relative to other heavy elements that is important because we 
interpolate the yield tables to the metallicity (i.e. mj,o/m t ot,o 
where the sum is over all elements heavier than helium) of the sim- 
ulation star particle. 

The ejected mass predicted by equation {7} can in principle 
be negative if the simulated relative abundance of some element is 
much lower than the relative abundance assumed in the yield calcu- 
lation and if the element is not produced in significant amounts. To 
prevent stellar particles from obtaining negative total abundances, 
we therefore impose that for each element j, the total ejected mass, 
which consists of the sum of the mass ejected by intermediate mass 
stars and SNe of types la and II, is non-negative at each time step. 
If it happens to be negative, we set it to zero. In practice we have, 
however, not observed this to occur in any of our simulations. 

If some w,j, e j has been set to zero, as in the scenario de- 
scribed in the previous paragraph, or if the simulated abundances 
(m Ji o:star/mtot,o : star) sim do not add up to unity (which can hap- 
pen when using smoothed abundances; see section l5"72l . then the 
sum of the ejected masses, ^ m j,ej, may differ from the total 
ejected mass, (m to t,cj)tbi, predicted by the yield tables. To ensure 
that we release the correct amount of mass, we therefore normalise 
the mass released from each element, including the total metal mass 
(see section |4T,l| >, such that the total mass released agrees with the 
value found in the yield tables. 



4.1.1 Tracking the total metal mass 

Of the many elements tracked by stellar evolution calculations, 
only a fraction contribute significantly to the radiative cooling 
rates. These elements tend to be the most abundant and thus the 
most observed. We therefore only track 9 elements individually: 
H, He, C, N, O, Ne, Mg, Si, and Fe (plus we take Ca and S, 
whose contributions to the cooling rates can be significant, see 
I Wiersma et all [2009, to be proportional to Si). Because we only 
track a finite number of elements, we cannot define the metal- 
licity to be Z = . nij/mtot where the sum is over all ele- 
ments heavier than He that are tracked by the simulation. Using 
Z = 1 — (mH + rriHe)/fTitot is also problematic since this defini- 
tion is susceptible to round-off errors. 

We therefore treat the total metal mass in the same manner 
as we treat individual elements, i.e., we include it as an addi- 
tional array, as has also b een done by othe r s (e.g. ILia et al l l2002l ; 
Kawata & Gibson] 120031: Ivaldarmnil [20031 ; [s ommer-Larsen et al.l 



20051 : 



Romeo et al 



2006). There are a number of ways to deter- 



produced. This is mainly due to non-conservation of nucleons and 
round-off errors (L. Portinari, private communication). We choose 
therefore to define Yz as the sum of the yields of all elements heav- 
ier than helium and redefine the hydrogen yield as: 
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mine the yield of the total metals released by a stellar population, 
Yz = rnz.p-.cy Since the sum of all Y, = mj tP - c j (including 
hydrogen and helium) should be zero (see eq. [6), one could sim- 
ply define Yz = — (Yh + Yu c ). Unfortunately, yield tables often 
give . Yj < 0, which implies that more mass is destroyed than 4.2 Results 



Figure 2. The fraction of the integrated mass ejected by an SSP that is due to 
SNII (solid), SNIa (dashed), and AGB (dot-dashed) as a function of the age 
of the SSP for two initial stellar metallicities: solar (black) and one per cent 
of solar (red). All curves assume a Chabrier IMF. The bottom-right panel 
corresponds to the total ejected mass, whereas the other panels correspond 
to individual elements, as indicated in the plot. While nearly all the ejecta 
come initially from SNII, for ages r > 10 s yr the contributions from AGB 
stars become significant for carbon, nitrogen and, if the metallicity is high, 
for other elements as well. The contribution from SNIa very important for 
iron for r > 10 9 yr. Note that the predictions are uncertain at the factor of 
two level, mainly due to the freedom in the normalization of the SNIa rate 
and the uncertainty in the yields. 



Note that Yz will differ from the sum of the metal yields of all 
elements tracked by the simulation if the yield tables contain more 
elements than are tracked by the simulation. The differences and 
round-off errors may not be significant in light of the uncertainties 
in our method (e.g., different yield sets, the definition of metallicity, 
etc.), but we choose to try to minimize the uncertainties at each 
step. 



Yn = -(Yn c + Yz) 



(8) 



Figures [2] and [3] show the result of combining all the ingredients 
described in appendix |A]using the implementation described in this 
section. Both figures provide information about the composition 
of the integrated (i.e., cumulative in time) ejecta of an SSP as a 
function of its age. While figure [2] shows the fraction of the mass 
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Figure 3. The abundances of various elements relative to iron in the integrated ejecta of an SSP as a function of its age for two metallicities: solar (solid) and 
one per cent of solar (dashed). All curves assume a Chabrier IMF. The abundances of all elements shown are initially enhanced relative to iron, particularly for 
lower stellar metallicities. The integrated abundances decrease relative to that of iron for r > 10 9 yr when large amounts of iron are ejected by SNIa. Note 
that the predictions are uncertain at the factor of two level, mainly due to the freedom in the normalization of the SNIa rate and the uncertainty in the yields. 



due to SNII, AGB and SNIa, figure[3]plots the abundance of a few 
elements relative to iron. 

At early times, r < 10 8 yr, the integrated ejecta are totally 
dominated by SNII and consequently the abundances relative to 
iron of elements such as carbon, nitrogen, oxygen, and silicon are 
highly supersolar. As the stellar population ages, the contributions 
of AGB stars become very important for carbon and nitrogen and 
those of SNIa for iron. Higher initial stellar metallicities tend to 
result in larger abundances relative to iron, although the effect be- 
comes smaller at late times because the yields of SNIa are indepen- 
dent of metallicity. 

Note that these figures show the properties of the mass re- 
leased integrated over time rather than at a particular time. For ages 
r > 10 8 yr all mass loss comes from AGB stars and SNIa and the 
elemental abundances relative to iron of the instantaneous ejecta 
are therefore much lower than those in the integrated ejecta that are 
shown here. 

We emphasize that the predictions shown in these figures are 
very uncertain, even for a fixed IMF. As we discuss in appendix |A3l 
the yields are uncertain at the factor of two level. For example, we 
rescaled the SNII yields of C, Mg, and Fe by factors of 0.5, 2, and 
0.5, respectively (note that these factors were included in Figs. [2] 
and[3] but not in any of the figures appearing in appendix [A3}. The 
relative contribution of SNIa is proportional to the normalization of 
the SNIa rate, which, as discussed in appendix IA4I is also highly 
uncertain. As shown in Fig.[2] this is particularly important for iron. 

We assumed a Chabrier IMF, but the results presented in this 
section would be similar for other IMFs as long as they have similar 
slopes for stellar masses above a few solar masses, as is for example 
the case for the Salpeter IMF. 



5 IMPLEMENTATION INTO SPH 
5.1 Enrichment scheme 

Having determined the ejected masses m-^ej for each stellar parti- 
cle during a given time step, we still need to define how much of 
this mass should go to each of the surrounding gas particles. We 
would like the mass transfer to be isotropic, but it is not obvious 
how to accomplish this. As discussed in detail in appendix A of 



IPawlik & Schavd < |2008|) , the transport of mas^j from a source par- 
ticle to its SPH neighbours is in general highly biased towards the 
direction of the centre of mass of the neighbours. In essence this is 
a consequence of the fact that particle-to-particle transport is only 
possi ble in directions where there are particles. 

IPawlik & Schavel (2008) demonstrate that usi ng an SPH ker- 
nel tra nsport scheme, as for example laid out by iMosconi et al.l 
results in a statistical bias that is much smaller than al- 
ternative schemes such as equal weighting or solid angle weight- 
ing. It should be noted, however, that while SPH weighting re- 
sults in only a small statistical bias, the transport from individual 
source particles is still highly anisotropic. As in previous works 
(e.g. I Steinmetz & Muellerl 19941: lKawatall200l|: iKawata & Gibson! 
2003 c iTornatore et al.ll2004l;|Scannapieco et alj|2005l : lstinson et al.l 
20061 : iKobavashi et alj|2007l : iTornatore et alj|2007h . we choose to 
employ SPH weights and transfer to each SPH neighbour fe of a 
given star particle the fraction of the ejected mass given by 



Wk 



(9) 



where h is the smoothing length of the star particle (which we de- 
termine in the same manner as for gas particles), r, is the distance 
from the star particle to neighbour i, W is the SPH kernel, and the 
su m in the denominator is over all SPH neighbours of the star parti- 
cle, [ibmatoreleral] ( 120071) have investigated changing the number 
of neighbours and kernel length used for distribution of heavy el- 
ements. They found that using more neighbours gives somewhat 
higher star formation rates, presumably because more particles are 
affected by metal-line cooling. We choose to use 48 neighbours, the 
same number as is used in the other SPH calculations. 

Finally, we note that we do not change the entropy or, equiv- 
alently, the internal energy per unit mass of the receiving gas par- 
ticles. The total internal energy of a receiving particle therefore in- 
creases when its mass increases. In reality the thermal energy of 
the gas surrounding a star undergoing mass loss will also change 
as a result of the mass transfer, but not necessarily by the amount 



6 IPawli k & Schaye 1 2008) considered photon transport, but the analysis is 
equivalent for mass transport. 
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that we implicitly assume. In fact, in the case of SN explosions we 
inject (part of) this energy explicitly (see Sj2}. On the other hand, 
many of the receiving gas particles may be star-forming, in which 
case we impose an effective equation of state (see ij2]l. Assuming all 
star particles lose 100 per cent of their mass and that all this mass is 
transferred to non-star-forming particles with T = 10 4 K, the total 
energy per unit stellar mass that we implicitly use to heat the trans- 
ferred mass to this temperature is only (c B /v w ) 2 /rj ~ 10~ 4 of the 
kinetic energy released by core collapse SNe, where c s is the sound 
speed of the ambient medium. Thus, the energy that we implicitly 
use in the mass transfer is negligible. 

Similarly, we also neglect the change in the momentum of star 
particles and their gaseous neighbours as a result of the change in 
the particle masses. We made this choice partly because we also do 
not attempt to follow the actual momentum of the ejecta except for 
core collapse SNe, in which case we inject the kinetic energy ex- 
plicitly. Neglecting momentum transfer from AGB stars and SNIa 
can be justified for the star particles if the mass transfer is isotropic, 
which it is in our case. For the receiving gas particles this argu- 
ment does not apply, even when averaged over many mass transfer 
events, if the stars and gas have systematically different bulk flows. 
However, even in that case the energy associated with this momen- 
tum transfer is negligible compared with the kinetic energy injected 
by core collapse SNe unless the bulk flow velocities of the stars dif- 
fer systematically from those of the gas by > 10 2 kms -1 . 

5.2 Smoothed metallicities 

The elemental abundances of a particle are not only of interest in 
the analysis of the simulations, they are also required during the 
simulation itself because gas cooling rates depend on them. Fur- 
thermore, the lifetimes and yields of star particles depend on their 
metallicities. 

Having tracked and transferred the total metal mass released 
by star particles during all time steps (see section |4.1.1| >, we could 
simply define the metallicity of a particle to be 

Z part EE ™. (10) 

m 

We will refer to metallicities computed using the above expression 
as "particle metallicities". Alternatively, we could define the metal- 
licity as the ratio of the SPH smoothed m etal mass density and 
the SPH smoothed g as mass density (as in lOkamoto et aT, I l2005l ; 
iTornatore et al.l2007l) . 

Z sm = ^L , (11) 
P 

where the gas and metal mass densities of particle i are given by 
the standard SPH expressions 

Pi = /J m 3 W{ri - rj, hi), (12) 
3 

pz,i = /^^mzjWjri -rj,hj), (13) 
3 

where the sum is over all SPH neighbours j and hi is the smoothing 
length of particle i. We will refer to metallicities computed using 
i ll It as "smoothed metallicities". Similarly, we will make use of 
particle and smoothed abundances of individual elements. Note that 
while particle abundances can only change when a particle gains or 
loses mass, smoothed abundances will generally vary from time 
step to time step because they depend on the distances between 
particles and their SPH neighbours. Although most studies do not 



explicitly define what they mean by "metallicity", particle metal- 
licities are generally used. 

The use of smoothed abundances has the advantage that it is 
most consistent with the SPH formalism. For example, in the ab- 
sence of metals gas cooling rates depend on p and it therefore 
makes sense for the metallic cooling rates to depend on pz ss 
Z STa p, computed in the same manner as p. Another advantage of 
smoothed abundances is that they partially counter the lack of metal 
mixing that is inherent to SPH. In particular, there will be many less 
particles with zero metallicity if smoothed abundances are used be- 
cause every neighbour of a particle with a non-zero particle metal- 
licity will have a non-zero smoothed metallicity. 

The reason SPH underestimates metal mixing is that, in the 
absence of some implementation of diffusion, metals are stuck to 
particles. Even in the absence of microscopic diffusion processes, 
this will result in a lack of metal mixing due to a straightforward 
sampling problem. This sampling problem is most easily illustrated 
by imagining a shell of constant thickness swept up by a spherically 
symmetric explosion driven by a single stellar particle in a uni- 
form, initially primordial gas, just after it has completed the trans- 
fer of metals produced by massive stars to its SPH neighbours (see 
Fig. |4j. As the metal-rich bubble expands, the local metal density 
should decrease as 1/r 2 . However, because the metals are stuck 
to a fixed number of particles, the local metal density will instead 
fluctuate strongly within the swept-up shell (Fig. [4] middle panel). 
The situation may actually be even worse in practice because we 
are using kinetic feedback which implies that only a subset of the 
SPH neighbours of the star particle are kicked (Fig. [4] right panel). 
Note that this sampling problem does not become smaller when the 
resolution is increased. 

The metal sampling problem of SPH can be reduced by imple- 
menting diffusion (and tuning the diffusion coefficient). We have, 
however, decided not to do this because diffusion is a physical 
process whereas the sampling problem is in essence a numerical 
problem. Attempting to solve it with a poorly constrained physical 
mechanism may have undesired consequences. Depending on the 
severity of the sampling problem, it may for example require an 
unphysical amount of diffusion. The addition of diffusion may also 
lead to the introduction of new physical scales (through the choice 
of a particular value for the diffusion coefficient). 

We therefore consider the problem of (turbulent) diffusion 
separate from the sampling problem described here. While diffu- 
sion processes may w ell be important and have a signi ficant im- 
pact on our results (e.g. lde Avillez & Breitschwerdtl2007l) . sub-grid 
recipes are required to implement them as they relevant physical 
processes are unresolved in cosmological simulations. 

Note that a static gas will not suffer from the sampling prob- 
lem but diffusion will still alter its metallicity distribution. A possi- 
ble way out may therefore be to implement a recipe for (turbulent) 
diffusion that dep ends explicitly on the relative velocities of neigh- 
bouring particles dGreif et al.l2008t) . Although this will still implic- 
itly use physical processes (e.g. hydrodynamical instabilities) to get 
around a numerical problem, it does have the important advantage 
that it asymptotes to the standard, no mixing, solution for the case 
of a static gas. We intend to test such sub-grid mixing prescriptions 
in the future. For now, we note that our use of smoothed abundances 
decreases the severity of - but does not solve - the sampling aspect 
of the mixing problem, leaving the physical aspect of the problem 
unaddressed. 

We decided to use smoothed abundances in the calculation of 
the cooling rates, the stellar lifetimes and the yields. When a gas 
particle is converted into a stellar particle its smoothed abundances 
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Figure 4. The enrichment sampling problem. A: A star particle enriches its neighbouring gas particles (red). B: The energy released by massive stars within the 
star particle drives its neighbours away. Because metals are stuck to particle the local metallicity in the shell fluctuates. C: Using kinetic feedback the problem 
is worse because only a small fraction of the neighbours are kicked. 



are frozen. Note that because the smoothed abundances of gas par- 
ticles are recalculated at every time step from the metal mass frac- 
tions, we need to store both smoothed and particle abundances and 
can thus compare the two. 

It is important to note that metal mass is only approximately 
conserved when smoothed abundances are used using the "gather 
approach", equation J13t . The non-conservation is exacerbated by 
the fact that the smoothed abundances of a particle are frozen when 
it is converted from a gas into a star particle. Metal conservation 
could be enforced by using a "scatter" approach, 

pz,i = ^ mzjW'jn - rj,hj), (14) 

3 

normalized such that J^. W" '{n — rj, hj) = 1 for every particle 
j. However, this is computationally more expensive and would in 
any case only conserve metal mass if the time steps of all particles 
were synchronized, which is not the case for GADGET. 

It is a well known problem that the SPH technique overes- 
timates densities of low density gas in regions wher e the density 
gradie nt is steep, which may result in overcooling iPearce et al.l 
dl999h . By using smoothed metal densities we are vulnerable to 
a similar problem that may result in an overestimate of the metal 
mass and hence increased overcooling. Indeed, in our low resolu- 
tion L025N128 and L100N128 runs the final total smoothed metal 
masses are, respectively, 3.6% and 18% greater than the final to- 
tal particle metal masses. However, as expected, the difference 
dereases rapidly with increasing resolution. For the L100N256, 
L100N512, and L025N512 runs the total smoothed metal mass is 
8.3% higher, 3.7% higher, and 0.14% lower than the total particle 
metal mass in the respective simulations. Thus, metal mass conser- 
vation errors are negligible for our highest resolution simulations. 

However, note that we compared total smoothed and parti- 
cle metal masses in simulations that always used smoothed abun- 
dances to compute the cooling rates, stellar lifetimes and yields. 
Because increased metal abundances will result in increased cool- 
ing and thus increased star formation and increased metal pro- 
duction, we expect the difference to increase if we compare the 
total smoothed metal mass predicted by a simulation employing 
smoothed abundances with the total particle metal mass predicted 
by a simulation employing particle abundances. We have tested 
this by comparing two L100N256 (and two L100N128) simula- 
tions. One simulation used smoothed abundances while the other 
used particle abundances for the calculation of the cooling rates. 
Indeed, we found that at z = the total smoothed metal mass in 
the L100N256 (L100N128) simulation using smoothed abundances 
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Figure 5. Stellar mass density as a function of redshift for 256 3 , 
100 h^ 1 Mpc runs. The black curve shows the total stellar mass density 
in L100N256 which employs smoothed abundances in the calculation of 
the cooling rates. The red curve indicates the total stellar mass in a simula- 
tion that used particle abundances for the cooling rates. The use of smoothed 
abundances increases the stellar mass produced by z = by 5 1 %. The blue 
curve used smoothed abundances for the cooling after multiplying them by 
a factor 1/1.08 to correct for the increase in the total metal mass, at fixed 
stellar mass, due to the use of smoothed abundances. The fact that the blue 
curve is close to the black one indicates that the increase in the stellar mass 
when using smoothed abundances is mostly due to the increased metal mix- 
ing rather than the small increase in the total metal mass. 



was 49% (82%) higher than the total particle metal mass in the cor- 
responding simulation using particle abundances. 

As can be seen from Figure [5] the simulation that uses 
smoothed abundances (black curve) produces about 1.5 times as 
many stars as the simulation that used particle abundances (red 
curve) for the calculation of the cooling rates. We performed a 
further test to determine if the increase in the stellar mass re- 
sulting from our use of smoothed abundances is due to the non- 
conservation of metal mass or to the increased metal mixing (i.e., 
metal cooling affects more particles). To this end we ran another 
version of the L100N256 simulation that used smoothed abun- 
dances but in which the smoothed abundances were multiplied by a 
factor 1/1.08, which is the ratio between the final total particle and 
smoothed metal mass in the simulation that used smoothed abun- 
dances for the cooling, before passing them to the cooling routine. 
As can be seen from Figure[5](blue curve), the result is much closer 
to the simulation using smoothed abundances (but without the re- 
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Figure 6. Two-dimensional probability density function of particle num- 
ber density in the smoothed metallicity - particle metallicity plane for the 
L100N512 simulation at z = 0. The contours are spaced logarithmically 
by 0.5 dex. The two metallicity definitions agree in highly enriched regions 
where the metals are well mixed, but at lower metallicities the smoothed 
metallicities are typically higher. The lowest possible smoothed metallic- 
ity for a particle with non-zero particle metallicity corresponds to the case 
when all its neighbours are metal-free and agrees with the lower envelope 
visible in the plot. 



duction factor) than to the simulation using particle abundances. 
Hence, the increase in the total stellar mass is mostly due to the 
increase in the metal mixing. 

This example demonstrates that the poor performance of SPH 
with regards to metal sampling and mixing is an important problem 
that can have a very large effect on the predicted star formation his- 
tory. On the other hand, we expect the differences to be smaller for 
our higher resolution simulations. The fact that non-conservation 
of metal mass is much less important than metal mixing and that 
the mismatch between the total metal masses becomes very small 
for our highest resolution simulations, support our decision to use 
smoothed abundances during the simulations. In the rest of this pa- 
per we will therefore only make use of simulations that employed 
smoothed abundances for the calculation of the cooling rates. 

In figure[6]we directly compare the 2 = smoothed and par- 
ticle metallicities in the L100N512 simulation. The contours show 
the logarithm of the particle number density in the Z sm — Z part 
plane. The dashed line indicates 1-1 correspondence. The two 
metallicities are strongly correlated for high metallicities, although 
significant scatter exists. This is expected, because high-metallicity 
gas will have been enriched by many star particles during many 
time steps and we therefore expect the metals to be well mixed. For 
low particle metallicities, however, the two metallicity definitions 
become essentially uncorrelated and the smoothed metallicities are 
typically higher than the particle metallicities. This is because a sin- 
gle enriched particle can give its neighbours a range of smoothed 
abundances depending on their distances. 

For a given particle metallicity, there is a well defined lower 
limit to the smoothed metallicity. This minimum smoothed metal- 



licity corresponds to that of a metal-bearing particle that is sur- 
rounded by metal-free particles. It is visible as the lower envelope 
of the contours parallel to the dashed line in figure [6] From equa- 
tions dl lb and J 1 3b we can see that the minimum smoothed metal- 
licity is Z sln .i — mz,iW(0, hi) I pi. Using the fact that for GAD- 
GET the value of the kernel at zero lag is W(0,hi) = 8/{nhf) and 
that for particles with rrii ~ m the kernel satisfies 47rpi/i 3 /3 ~ 
iVngbm (where the number of SPH neighbours -/V ng b = 48 in our 
simulations), it is easy to show that Z sm .i ~ 2Z part ,i/9. 

In summary, we use smoothed abundances during the sim- 
ulation because it is consistent with the SPH method and be- 
cause it reduces, though not eliminates, the metal sampling prob- 
lem. Smoothed and particle metallicities are very close in high- 
metallicity regions but in low-metallicity regions the smoothed 
abundances will typically exceed the particle ones. Using SPH sim- 
ulations to study low-metallicity gas and stars is hence problematic 
because the results are sensitive to the definition of metallicity that 
is used. The use of smoothed abundances increases the number of 
particles for which metal cooling is important and thus the star for- 
mation rate, particularly when the resolution is low. 



6 THE PREDICTED DISTRIBUTION OF METALS 

In this section we will investigate the cosmic metal distribution pre- 
dicted by our simulations. We will only present results from the 
simulations listed in Table [2] which are drawn from the reference 
simulations of the OWLS set. While the simulations analysed here 
differ in terms of their box sizes and particle numbers, they all used 
the same simulation code and the same prescription for star forma- 
tion, galactic winds and cooling. In a future publication we will use 
the full OWLS set to study the effect of changes in the baryonic 
processes on the metal distribution, which turns out to be very im- 
portant. Here we will instead focus on convergence tests and on the 
effect of different metallicity definitions. 

There are a number of approaches we can take to study the 
cosmic metal distribution in our simulations. Before analysing the 
metal distribution at z — in more detail, we will investigate the 
evolution of the metal mass fractions in various components. We 
will then proceed to discuss the dependence of the results on the 
definition of metallicity ( ij6.lt . the size of the simulation box (ap- 
pendix|Bj and the numerical resolution (appendix[Cj- While appen- 
dices [B] and [C] will make use of the full suite of simulations listed 
in Table [2] we will only consider the L100N512 (which contains 
2 x 512 3 particles in a box that is 100 hT 1 Mpc on a side) run 
in the rest of this section. We note, however, that our much higher 
resolution L025N512 simulation (which was, however, stopped at 
z — 1) gives qualitatively similar results. 

Figure[7]shows the evolution of the fraction of the metal mass 
contained in various components. By isolating the gaseous met- 
als into various phases, we can develop a picture of where the 
metals are at a given redshift. To that end we have plotted the 
total fraction of the metal mass in stars (left, solid), star-forming 
gas (left, dashed), and non-star-forming gas (left, dot-dashed). Re- 
call (see [j2]( that by star-forming we mean gas particles that are 
allowed to form stars and for which we impose an equation of 
state. These particles all have densities «h > 10 _1 cm~ 3 and 
can be thought of as ISM gas. We have divided the non-star- 
forming gas into cold-warm (T < 10 J K; middle) and warm-hot 
(T > 10 5 K; right) components. The cold-warm gas has been fur- 
ther subdivided into diffuse (p < 10 2 (p); middle, solid) and halo 
(p > 10 2 (p); middle, dashed) gas and we have subdivided the 
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Figure 7. Evolution of fractional metal mass in various components for the L100N512 simulation. Left: Stars (solid), star-forming gas (dashed) and non- 
star-forming gas (dot-dashed). Note that these curves add up to unity. Centre: Cold-warm IGM (non-star-forming, p < 100{p), T < 10 5 K) (solid), and 
cold-warm halo gas (non-star-forming, p > 100{p), T < 10 5 K) (dashed). Right: WHIM (non-star-forming, T > 10 5 K) (solid) and ICM (non-star-forming, 
T > 10 7 K) (dashed). The parts of the curves corresponding to redshifts for which the total metal mass is smaller than 10~ 6 of the total baryonic mass have 
been omitted because they become very noisy. While most metals initially reside in star-forming gas, by z = most are locked up in stars. At present most of 
the gaseous metals reside in the WHIM. 




Figure 8. Evolution of the mass-weighted metallicities of various components for the L100N512 simulation. The panels and line styles are identical to 
those shown in Fig. [7j The metallicity of stars and the ISM evolves much more weakly than that of cold-warm diff use gas. The mean m etallicity of the 
WHIM is ~ IP -1 Z p) at all times. Symbols indicat e observational estimat es of metallicities in vari ous phases - stars: jpallazzi et a l. 2008) (pl us sign) and 
lHallidav et al. 2008) (c ross symbol); col d halo gas: JProchaska et al.ll2003l) (squares); diffuse IGM: jAguirre et al]|2008l) (triangle) and Ischave et al 2003) 
(diamond); ICM: ISimionescu et al.l l2009) (orange region). 



warm-hot, non-star-forming gas into the warm-hot IGM (WHIM) 
(10 5 K < T < 10 7 K; right, solid), and the intracluster medium 
(ICM) (T > 10 7 K; right, dashed). 

Focusing first on the left-hand panel of Fig. [7] we see that the 
most prominent trend is the strong decrease in the fraction of the 
metals in star-forming gas. While this phase harbors most of the 
metals at z = 8, by z = it contains only a small fraction. The 
opposite is true for the stars: while most of the metals are initially 
in the gas phase, today the majority of the metals are locked up in 
stars. The metal fraction contained in non-star-forming gas varies 
between twenty and forty per cent. 

iDave & Oppenheimeil j2007h have recently used a cosmologi- 
cal SPH simulation (using 256 3 gas particles in a 32 h' 1 Mpc box) 
to investigate the cosmic distribution of metals. While we postpone 
a detailed comparison to a future pap er, we note a clear, qualita- 
tive difference in the results. While Pave & Op penheimed d2007l) 
find that initially (at z = 6) the diffuse IGM contains much more 
metal than the ISM, we find the opposite. We will show elsewhere 
that this difference mostly refl ects the higher initial mass lo ading 
of galactic winds assumed by Dave & Oppenheimer (2007). This 
difference illustrates the importance of varying the subgrid recipes 
used in the simulations. We will report on such variations else- 
where. 

The other two panels of Fig.|7]show that while the metal frac- 
tion in cold-warm gas decreases with time, the fraction contained in 
warm-hot and hot gas increases. This is not unexpected, as growth 
of structure will result in more and stronger accretion shocks, heat- 
ing an increasing fraction of the gas to high temperatures. While 



most of the metals in non-star-forming gas are initially cold, more 
than half of them reside in the WHIM for z < 2. The ICM, on the 
other hand, still contains only a few percent of the metal mass by 
z = 0, an order of magnitude less than the WHIM. 

Besides the metal mass fractions, the metallicities of the vari- 
ous components are also of interest. Together with the evolution of 
the baryonic mass fractions, which we will present elsewhere, the 
metallicities determine the metal mass fractions. Figure [8] shows 
the evolution of the mass-weighted mean metallicities of the same 
components that were shown in Fig. [7] 

As expected, the metallicity of the star-forming gas roughly 
traces the metallicity of the stars throughout the simulation, al- 
though the ratio of the mean stellar metallicity to the mean ISM 
metallicity decreases slightly with time. The stellar and ISM metal- 
licities increase relatively slowly with time. At z — 3 the mean 
stellar metallicity is about 0.5 dex lower than at z = 0. Interest- 
ingly, the WHIM and the ICM have a metallicity ~ 10~ 1 Z Q at all 
times. This lack of evolution is in stark contrast to the behavior of 
the cold-warm, diffuse gas whose metallicity increases rapidly with 
time. Consequently, the differences in the metallicities of the vari- 
ous components become smaller at lower redshifts. However, even 
by z = the metallicity of the diffuse, cold-warm IGM is only of 
order 1 per cent of solar. 

Although our pimary aim is not to perform a direct compar- 
ison to observations, we have included in Fig. [8] a number of ob- 



© 0000 RAS, MNRAS 000, 000-000 



12 R. P. C. Wiersma et al. 




Figure 9. Gas distribution weighted by mass (black) and volume (red) in the temperature - density (left), metallicity - density (centre), and metallicity - 
temperature (right) planes at z = for the L100N512 simulation. Note that star-forming gas has been excluded from the right-hand panel. The contours 
are spaced by 1 dex. While low-density gas displays a wide range of metallicities, there is a strong positive gradient of metallicity with density that becomes 
narrower for higher densities. Combined with the fact that low-density gas must account for a larger fraction of the volume than the mass, the metallicity- 
density correlation implies that volume-weighting favors lower metallicities than mass weighting. There is no well defined relation between temperature and 
metallicity. 
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Figure 10. Gas mass (black) and gas metal mass (red) distribution in the temperature - density (left), metallicity - density (centre), and metallicity - temperature 
(right) planes at z = for the L100N512 simulation. Note that star-forming gas has been excluded from the right-hand panel. The contours are spaced by 
1 dex. While gaseous metals track the general gas mass distribution at high densities (or metallicities), at the densities corresponding to diffuse structures 
(p < 10 2 (p)) metals reside in gas that is both hotter and more metal-rich than is typical for these densities. 



servational estimates of metallicitjQ in various phases. In the left 
panel we show the esti mates for the mean s tellar metallicity at 
z = and z = 2 from iGallazzi et all (200$ and lHallidavetal.1 
(2008), respectively. Our simulations fall about 0.2 dex above the 
redshift zero point and land abou t 0.3 dex abov e the h igher red- 
shift point. We note however, that lHallidav et al] 1 20081) find their 
result to be low compared to that of Erb et alf i feOOrj) They ex- 
plaine d that this could represent an alpha enha ncement at high red- 
shift jErb et alj|200r3 measured oxygen while lHallidav et al ||2008| 
measured iron). We compare our predictions for overdense, cold, 
non-sta rforming gas to the obse rvations of damped Lya (DLA) sys- 
tems of Procha ska et alj j2003h and find excellent agreement. Our 
measurement of the mean metallicity of the z ~ 3 IGM agrees 
with the measurement based on carbon of lSchave et alj d2003l), bu t 
falls slightly below that based on oxygen of lAguirre et alj (2008). 
In appendix [C] we show that the metallicity of this phase is partic- 
ularly sensitive to resolution and that higher resolution simulations 
may predict higher diffuse IGM metallicities. Finally, the metallic- 
ity of the low-redshift ICM falls along the lower ran ge of cluster 
measurements compiled by ISimionescu et aD (|2009). For this we 
compiled a range of all metallicity measurements in the most outer 
regions of the clusters (regions which should dominate the mass av- 
eraged metallicity). We note that since most radial profiles flatten 



The observations were converted to our solar abundances where neces- 
sary. 



at outer radii, this should be a safe assumption. From these compar- 
isons we conclude that our simulation is in reasonable agreement 
with a range of observations. 

Before investigating the robustness of these results to numer- 
ical effects, we will study the z = metal distribution in more 
detail. To interpret the results, it is, however, helpful to first con- 
sider the gas mass rather than the metal mass distribution. 

The solid contours in the left-hand panel of figure [9] indicate 
the z = gas mass-weighted 2-dimensional (2-D) probability den- 
sity function (PDF) in the temperature-density plane. That is, the 
contour values correspond to dM s / (AI g ,totd log -A-d log T). Sim- 
ilarly, the red contours show the volume-weighteqj 2-D PDF, i.e., 
dV/(Vtotdlog -A-dlogT). All contour plots in this section show 
contours spaced by one dex, with each plot containing the same 
levels for both linestyles. 

There are a number of distinct features in this temperature- 
density diagram. Much of the mass, and nearly all of the volume, 
resides in a narrow strip ranging from about (log p/ (p) , log T) ~ 
(—1.5, 3) to ~ (1.5, 4.5), which corresponds to the diffuse, photo- 
ionized IGM. The relation between temperature and density of this 
component is set by the balance between photo-heating and adia- 



8 Volume weighting was implemented as weighting by h f / ^ . h^, where 
hi is the SPH smoothing length of particle i. We verified that using m/p 
instead of h? gives nearly identical results. 
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batic cooling faui&Gnedinl 19971) . For gas at higher overdensities 
photo-heating is balanced by radiative rather than adiabatic cooling 
and the temperature slowly decreases with increasing density. The 
plume of much hotter gas corresponds to gas that has been shock- 
heated, either through gravitationally induced accretion shocks or 
galactic winds. For gas with densities pj (p) > 10 there is less 
mass at temperatures T ~ 10 K than there is at lower and higher 
temperatures because the radiative cooling rates peak at these tem- 
peratures (due mostly to collisional excitation of oxygen, carbon, 
and helium), resulting in a bimodal distribution of the temperature 
in that density regime. The well-defined T-p relation at the highest 
densities reflects the equation of state P oc p^J z ', which we im- 
pose on gas with densities exceeding our star formation threshold 
of Tin =0.1 cm" 3 (pj (p) ~ 10 6 at z = 0). As discussed in iJU 
the temperature of this gas merely reflects the imposed effective 
pressure of the unresolved multiphase ISM. 

The middle and right panels of Fig.|9]show how the gas mass 
(black) and the volume (red) are distributed in the metallicity- 
density and the metallicity-temperature planes, respectively. We 
have excluded star-forming gas from the right-hand panel since 
its temperature would merely reflect the pressure of the equation 
of state that we have imposed. We will exclude this component 
from all plots investigating temperature dependencies, except for 
temperature-density diagrams such as the left-hand panel of Fig. [9] 
(because in that case the gas density can be used to identify the 
star-forming gas). 

While there is a clear positive correlation between metallic- 
ity and density for highly overdense gas (see the black contours in 
the middle panel), gas at low overdensities can have a large range 
of metallicities extending all the way to zero (outside the plotted 
range). Note that the highest density gas is highly supersolar. This 
reflects the fact that massive galaxies in this simulation have too 
high metallicities because galactic winds cannot escape their high- 
pressure ISM. In a future paper we will show that this problem 
is not present for some other implementations of galactic winds 
driven by SNe and, particularly, for more efficient feedback mech- 
anisms such as energy injected by active galactic nuclei. The right- 
hand panel shows that there is no well-defined relation between 
metallicity and temperature. Gas at a given temperature displays a 
wide range of metallicities. 

Figure [Tol uses the same panels as Fig. [9] to compare the dis- 
tribution of gas mass (black) and metal mass (red). Since we use 
smoothed metallicities, metal mass is defined as Z sm M & . Note that 
the black contours are identical to those shown in Fig. [9] Focusing 
first on the left-hand panel, which shows the distributions in the 
temperature-density plane, we see that while metal and gas mass 
track each other well for high overdensities (although the gas mass 
is distributed over a narrower range of temperatures for a given 
density), they differ in one important respect for pj (p) < 10 2 . 
Whereas most of the mass is concentrated at the lowest densities 
and temperatures, this diffuse, photo-ionized IGM contains only a 
small fraction of the metals. Instead, the vast majority of metals 
in low-density gas are hot (T > 10 K). This is probably caused 
by the fact that high- velocity winds are required to transport met- 
als to regions far from gala xies and that such w inds shock-heat the 
mediu m. This agrees with iTheuns et alj J2002h and lAguirre et al.l 
(2005), who also found that most of the diffuse metals reside in 
hot gas, using different implementations of feedback from star for- 
mati on and using s imulations that did not include metal-line cool- 
ing. lAguirre et alj d2005h speculated that the inclusion of radia- 
tive cooling by metals might be important, but our simulations 
show that even when metal-line cooling is included, the WHIM 



(5 < log T < 7, pi (p) < 10 3 ) is an important reservoir of cos- 
mic metals. From the right-hand panel of Fig. [9] we can see that 
although some WHIM gas has very low metallicities, much of it 
has Z > 10" 1 Z Q . 

The left-hand panel of Fig. QT] shows the cumulative 
PDF for the gas metal mass as a function of logp/ (p), i.e., 
T7 — f, P dp' d ( M s Zs ™'> w here Mz„,„ t „ t is the total smoothed 

Af Zam>t ot JO r dp' "sm.tot 

metal mass in gas. Similarly, the second and third panels from the 
left show the cumulative PDF of the gas metal mass as a function of 
temperature and metallicity and the right-hand panel shows the cu- 
mulative PDF of the stellar metal mass as a function of metallicity. 
While the gas metal mass is spread over remarkably large ranges 
of densities and temperatures, it is concentrated in relatively high- 
metallicity gas. Note, however, that metallicity is evaluated locally, 
at the resolution limit of the simulation. High metallicity gas there- 
fore includes patches of enriched gas embedded in structures whose 
overall metallicity is low. 

Considering one variable at a time, ninety percent of the 2 = 
gas metal mass resides in gas with 10 -0 - 5 < pi (p) < 10 7 - 5 , 
10 4 K < T < 10 65 K, or KT 1 - 5 < Z am /Z® < 10 07 . Half 
the gaseous metal mass resides in gas with p < 10 2 (p) which im- 
plies that diffuse and collapsed structures contain similar amounts 
of metal mass. In terms of temperature the midpoint lies at about 
10 5 K, which means that half the gas metal mass resides in shock- 
heated gas. Clearly, any census aiming to account for most of the 
metal mass in gas will have to take a wide variety of objects and 
structures into account. 

In the next section we will investigate the metal distribution 
and its evolution in more detail, including the sensitivity of the re- 
sults to the definition of metallicity. The convergence of the predic- 
tions with respect to the size of the simulation box and the numeri- 
cal resolution is studied in appendices [B] and [C] respectively. 



6.1 Smoothed vs. particle metallicity 

As discussed in i]5.2l there is some ambiguity in the definition of 
the metallicity in SPH. We have opted to use SPH smoothed metal- 
licities Z sm = pz/p rather than particle metallicities 
mz/m, because using smoothed abundances is consistent with the 
SPH method and because it partially counters the lack of metal mix- 
ing that results from the fact that metals are stuck to particles. We 
already demonstrated in i]5.2l that the definition of metallicity is 
particularly important for low metallicity gas and that it has impor- 
tant consequences for gas cooling rates (and hence for the predicted 
star formation histories). Here we investigate its effect on the metal 
distribution by comparing the gas metal mass distributions com- 
puted using smoothed and particle metallicities. Note, however, that 
smoothed metallicities were always used during the simulation for 
the calculation of the radiative cooling rates, stellar lifetimes, and 
stellar yields. 

Figure [12] compares the gas metal mass distribution 
using smoothed (black) and particle (red) metallicities in 
the temperature-density (left), metallicity-density (middle) and 
metallicity-temperature (right) planes. Note that the black contours 
are identical to the red contours in Fig.[l0] The smoothed and parti- 
cle metal distributions trace each other reasonably well, except for 
the diffuse, photo-ionized IGM (p/ (p) < 10 2 , T < 10 5 K; left 
panel) and generally for low-metallicity gas (Z <C Zq ; middle and 
right panels). When smoothed metallicities are used larger fractions 
of the metal mass reside in low-metallicity gas, which is consistent 
with Fig. [6] 
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Figure 11. Cumulative probability density function of the gas metal mass as a function of gas density (left), temperature (centre-left), and metallicity (centre- 
right) at z = for the L100N512 simulation. The vertical, dotted line in the left panel indicates our star formation threshold. Also note that star-forming gas 
has been excluded from the temperature plot. The right-hand panel shows the cumulative PDF of the stellar metal mass as a function of the stellar metallicity. 
While the gas metal mass is distributed over a wide range of densities and temperatures, it is concentrated in gas that has a relatively high (local) metallicity. 




Figure 12. Gas smoothed metal mass (black) and gas particle metal mass (red) distribution in the temperature - density (left), metallicity - density (centre), 
and metallicity - temperature (right) planes at z = for the L100N512 simulation. Note that star-forming gas has been excluded from the right-hand panel. 
The contours are spaced by 1 dex. Except for low-metallicity gas, which mostly has low densities and temperatures, smoothed and particle metallicities yield 
similar results. 



The left-hand panel of figure [O] shows the 1-D 
gaseous metal mass weighted PDF for the gas density, i.e., 
d(M g Z)/ (Mz, totdlog as a function of logp/{p), where 
Z is the smoothed/particle metallicity and Mz.tot is the total 
smoothed/particle metal mass in gas. Thus, if the y-axis were 
linear (which it is not), the fraction of the area underneath the 
histogram that is in a given bin would correspond to the fraction 
of the gas metal mass in that bin. Similarly, the second and third 
panels from the left show the metal-mass weighted PDF of the 
temperature and metallicity, respectively, and the right-hand panel 
shows the metal mass-weighted PDF of the stellar metallicity. The 
distribution of metal mass over density is similar for smoothed and 
particle metallicities, but using smoothed metallicities results in a 
significant shift of metal mass from T ~ 10 6 K to T < 10 5 K. 
While the differences in the metal distributions as a function 
of metallicity are small for stars, presumably because they tend 
to form in well-mixed, high density gas, the differences are 
very substantial for the gas. As expected, the use of smoothed 
metallicities results in a prominent shift of gaseous metal mass to 
lower metallicities. For stars, the particle metallicity PDF exceeds 
that for smoothed abundances by about an order of magnitude in 
the range —4 < log Z < — 2 (i.e., outside the plotted range). 

We stress, however, that because the metallicities are evalu- 
ated locally, i.e. at the spatial resolution limit, they may differ sub- 
stantially from the overall mean metallicity of the structure con- 
taining the particles. For structures that contain many particles the 
mean metallicity will be insensitive to the definition of metallicity 
that is used. On the other hand, local gas metallicities are the ones 



that are used for the calculation of the cooling rates, so they are 
certainly important. 

In summary, SPH simulations such as the ones presented here 
are not well-suited to study the metal distribution in low-metallicity 
gas. Because metals are stuck to particles, SPH suffers from a sam- 
pling problem which becomes worse in regions where only a small 
fraction of the particles have been enriched. Consequently, in the 
low-metallicity regime the distribution of mass is sensitive to the 
definition of metallicity. However, the bulk of the metals reside in 
gas of relatively high metallicity for which the predicted metal dis- 
tribution is not strongly affected by the lack of metal mixing. Ex- 
cept for the small fraction of stars that have Z < 10 _1 Z@, the 
results are more robust for stars because they form in high-density 
gas that is better mixed. This is probably because high-density gas 
particles tend to be near star particles which means that they are 
likely to have received metals during many time steps. 



7 SUMMARY 

The elemental abundances of stars, galaxies, and diffuse gas are 
of great interest because they determine radiative cooling rates and 
the stellar initial mass function, because they can give insight into 
physical processes such as the interactions between galaxies and 
the IGM, because they change the strength of lines that are used 
as diagnostics of physical conditions and because they determine 
the observability of lines used as tracers of various gas phases. In 
this paper we have presented a method to follow the timed release 
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Figure 13. Probability density function of the gas density (left), temperature (centre-left) and smoothed/particle metallicity (centre-right) weighted by 
smoothed (black) and particle (red) metal mass. The black (red) histogram in the right-hand panel shows the PDF of the smoothed (particle) stellar metallicity 
weighted by smoothed (particle) metal mass. All panels correspond to the L100N512 simulation at z = 0. The vertical, dotted line in the left panel indicates 
the threshold for star formation. Note that star-forming gas was excluded from the temperature panel. The predictions for smoothed and particle metallicities 
only differ significantly for low-metallicity gas. 



of individual elements by stars, including mass loss through stellar 
winds and SNe, and their subsequent dispersal through space. 

The ingredients of our stellar evolution module include a 
choice of IMF (we use Chabrier), stellar lifetimes as a function 
of metallicity (which we take from IPortinari et aTll 19981) . the rate 
of type la SNe for an SSP as a function of time (we use an em- 
pirical time delay distribution tuned to obtain agreement between 
the evolution of the observed, cos mic rates of SN Ia and star for- 
mation), and stellar yields (we use|Marigo 2001 for AGB stars, 
Portinari et al. 1998 for core collapse SNe, and the W7 model of 
Thielemann et : alj2003l for SNIa). 

We compared different sets of nucleosynthetic yields taken 
from the literature after integrating them over the IMF. While the 
different studies predict similar total ejected masses, the predictions 
for heavy elements differ substantially. In particular, the ejected 
masses of individual elements differ by factors of a few for AGB 
stars and - for elements heavier than nitrogen - for core collapse 
SNe. Thus, even disregarding the fact that most nucleosynthetic 
yield calculations ignore potentially important effects such as ro- 
tation, the elemental ratios are uncertain by factors of a few for a 
fixed IMF. Predictions for abundances relative to iron are even less 
robust because much, perhaps even most, of the iron released by an 
SSP is produced by SNIa and the normalization of the SNIa rate is 
only known to within factors of a few. 

Our implementation of mass transfer implicitly splits, for each 
element, the ejected mass into terms accounting for the mass that 
is produced (minus destroyed) and the mass that is simply passing 
through. The latter term is assumed to scale with the initial abun- 
dance, allowing us to correct for relative abundance variations with 
respect to solar at a fixed metallicity0 We have taken special care 
to correctly track the total metal mass even if not all elements are 
tracked explicitly. 

We discussed two possible definitions of metallicity of SPH 
particles: the commonly used "particle metallicity", defined as the 
ratio of the metal to the total mass (mz/m), and the "smoothed 
metallicity", defined as the ratio of the metal mass density to the 
total mass density (pz/p), where each density is computed using 
the standard SPH formalism. We argued that smoothed metallicities 
(and, analogously, smoothed elemental abundances) are preferable 
because they are most consistent with the SPH formalism (particu- 



9 Note that this treatment still ignores the fact that the nucleosynthesis may 
itself depend on the relative abundances. 



larly with regards to radiative cooling rates which depend explicitly 
on densities) and because they partially counter the lack of metal 
mixing that is inherent to SPH. We discussed in some detail the fact 
that SPH underestimates metal mixing because metals are stuck to 
particles and how this results in a sampling problem that may not 
become smaller when the resolution is increased. While the use 
of smoothed abundances eases this problem somewhat, it by no 
means eliminates it entirely (while also leaving the potential for 
small scale metal mixing due to physical processes unaddressed). 

The use of smoothed abundances (which are frozen when a 
gas particle is converted into a star particle) leads to a slight non- 
conservation of metal mass, but we showed this to be negligible 
for our high-resolution simulations. A comparison of smoothed 
and particle metallicities (in a simulation that used smoothed abun- 
dances for the cooling rates and the stellar evolution) shows that 
while they are similar at high metallicities, they differ strongly at 
low metallicities, with smoothed abundances typically exceeding 
particle abundances by large factors. In particular, there are many 
more particles with non-zero smoothed metallicities than with non- 
zero particle metallicities. Since there is no overwhelming reason 
to prefer one definition of metallicity over the other, the fact that the 
choice matters indicates that care must be taken when interpreting 
the predictions of SPH simulations at low metallicities. 

A comparison between two low-resolution simulations, one 
using smoothed abundances and the other using particle abun- 
dances for the calculation of the cooling rates, revealed large differ- 
ences in the star formation rate. By redshift zero the simulation em- 
ploying smoothed abundances had formed about 1.5 times as many 
stars as the run using particle metallicities. We demonstrated that 
this difference could not be attributed to non-conservation of metal 
mass and was instead caused by the increased mixing (i.e., metal 
cooling is important for more particles when smoothed abundances 
are used). Although the difference is expected to decrease with in- 
creasing resolution, we conclude that it will be necessary to solve 
the metal mixing problem before SPH simulations can be used to 
make precise predictions for the cosmic star formation rate. 

We used a suite of large simulations with up to 2 x 512 3 par- 
ticles to investigate the distribution of heavy elements. All simula- 
tions used identical physical parameters and sub-grid modules. The 
simulati ons made use of recently developed modules for star for- 
mation dSchave & Dalla Vecchiall2008h and kinetic feedback from 
core collapse SNe jDalla Vecchia & Schavej 120081) . We followed 
all 11 elements (H He, C, N, O, Ne, Mg, Si, S, Ca, Fe) that 
IWiersma et all (2009) found to contribute significantly to the cool- 
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ing of photo-ionized plasmas. Radiative cooling (including photo- 
heating) was implemented element-by-element, assuming the gas 
to be optically thin and in ionization equilibrium in the pres- 
ence of an evolving photo-ionizing radiation background. We did 
this by making use of tables tha t had been pre-comp uted using 
CLOUDY following the methods o flWiersma et al. ( 2009). Element- 
by-element cooling and the inclusion of photo-ionization - not only 
for H and He, but also for heavy elements - are both novel features 
for cosmological simulations. 

Our predictions for the metallicities of various baryonic 
phases are in reasonably good agreement with available observa- 
tions. While most of the cosmic metal mass initially resides in gas 
of densities typical of the ISM (na > 10 _1 cm -3 ), by redshift 
zero stars have become the dominant reservoir of metal mass and 
the ISM contains only a small fraction of the metals. Diffuse gas 
(tih < 10" 1 cm -3 ) contains a significant fraction of the metals at 
all times. Except at very high redshifts, most of the diffuse metals 
reside in the WHIM (10 5 K < T < 10 7 K), but the metal mass re- 
siding in the ICM (T > 10 7 K) is always negligible. By the present 
time, the gaseous metal mass is distributed over a wide range of 
densities (5, 50, and 95 per cent of the metal mass resides in gas 
with pi (p) < 10 -0 ' 5 , 10 2 , and 10 75 , respectively) and tempera- 
tures (5, 50, and 95 per cent of the metal mass resides in gas with 
T (K) < 10 4 , 10 s , and 10 6 ' 5 , respectively), but is concentrated 
in gas that has a relatively high, local metallicity (95 per cent has 
Z > 10 -1,8 Z e averaged over the mass scale corresponding to the 
resolution limit). Clearly, any census aiming to account for most 
of the metal mass will have to take a wide variety of objects and 
structures into account. 

Although the mean stellar metallicity slowly increases with 
time, it is already of order ten per cent of solar by z = 8. The evolu- 
tion of the metallicity is much stronger for cold- warm (T < 10 5 K) 
diffuse gas, particularly for p/ (p) < 10 2 . By redshift zero the dif- 
fuse, photo-ionized IGM has a metallicity ~ 10~ 2 Z®, while the 
cold-warm gas in halos has Z > 10" 1 Zq . Interestingly, the metal- 
licity of the WHIM and the ICM is ~ 10" 1 Z Q at all redshifts. 

A comprehensive convergence study (see appendices [B] and 
[Cl revealed that, except for the ICM, a 50 ft -1 Mpc box is suf- 
ficiently large to obtain a converged result for the cosmic metal 
mass fractions and metallicities down to z = and 12.5 /i -1 Mpc 
suffices for z > 2. The convergence with respect to numerical res- 
olution was consistent with expectations based on a comparison 
with the Jeans scales. Our use of an effective equation of state for 
star-forming gas guarantees that the Jeans mass does not fall below 
/g 2 10 7 h -1 Mq, where / g is the local gas fraction, but even this 
scale is only marginally resolved in our highest resolution simula- 
tions (which have m g ~ 10 6 M®). We found that it is much 
easier to obtain converged predictions for the metallicity of the dif- 
ferent components than for their metal mass fractions. While the 
simulations presented here can provide (numerically) robust pre- 
dictions for the metallicities at z < 4, higher resolution models 
may be required for higher redshifts. 

Here we have investigated only a single set of physical param- 
eters. We postpone such a comparison to a follow-up paper that will 
investigate how the predictions for the metal distribution vary if we 
change the physical assumptions and sub-grid prescriptions, such 
as the cosmology, the star formation recipe, the implementation and 
efficiency of galactic winds, the cooling rates, and feedback from 
AGN. In this way we hope to isolate the processes that drive the 
evolution of the cosmic metal distribution. 
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APPENDIX A: UNCERTAINTIES AND DIVERSITY IN 
STELLAR EVOLUTION CHOICES 

In this appendix we explain in detail our implementation of stellar 
evolution, giving an indication of possible sources of uncertainties 
in chemodynamical simulations. 

We need a method for recycling products of stellar evolution 
back into the ISM. We therefore take the standard assumption that 
a star particle represents a simple stellar population. The feedback 
processes originating from stellar evolution that we will consider 
are winds from asymptotic giant branch (AGB) stars, type la super- 
novae (SNe), type II (i.e. core collapse) SNe and the winds from 
their progenitors. The above stellar processes differ in a number 
of ways, such as their nucleosynthetic production, their energetic 
output, and the timescale on which they act. 

Progenitors associated with AGB stars are typically intermedi- 
ate mass (0.8 M Q < M < 8 M Q ) stars wh ich have long (10 s yr < 
r < 1 10 yr) main sequence lifetimes (see lKippenhahn & Weigertl 
1 19941) . For AGB stars the mass loss occurs during a period that is 
very short compared with both their lifetimes and the dynamical 
timescales in cosmological simulations and we therefore assume 
the mass loss for an AGB star of given mass to happen within a 
single simulation time step at the end of the main sequence life- 
time. The kinetic energy injected by AGB stars is assumed to be 
unimportant, which is an excellent approximation both because ob- 
served AGB wind velocities are not large compared to the velocity 
dispersion in the ISM and because it takes an SSP billions of years 
for all the energy to be released. AGB stars are major producers of 
carbon and nitrogen, both of them being among the few elements 



that can be detected in the diffuse IGM (e.g. ISchave et aT]|2003l: 
iFechner & RichteJl2008l) . 

Since type II SNe (S NII) only result from massive stars (see 
ICrowther & Smart! 120071 and references therein), most of these 
events take place within tens of millions of years after the birth 
of a stellar population. For SNII progenitors mass loss on the main 
sequence can be substantial. However, since their main sequence 
lifetimes are short compared with the dynamical timescales in cos- 
mological simulations, it is still a good approximation to release all 
the mass ejected by stars of a fixed initial mass in a single time step 
at the end of their lifetimes. Note that if the time step governing a 
star particle is shorter than the lifetime of the least massive SNII 
progenitor, then the mass will still be ejected over multiple time 
steps. 

Type II SNe explosions dump so much energy in the ISM 
(~ 10 51 erg of kinetic energy per SN) in a short time that they 
may drive galactic- scale outflows from starbursting galaxies (e.g. 
IVeilleux et al.ll2005h . The implementation of these winds is trou- 
blesome since simulations lack the reso l ution to implement them 
corre ctly (e.g. ICeverino & K lypin 2007|; iDalla Vecchia & Schavd 
2008). In our simulations forty percent of the kinetic energy from 
SNII is injecte d in kinetic form as discussed in section [2] and in 
more detail in Dalla Vecchia & Schavel d2008l) . The remainder is 
assumed to be lost radiatively on scales below the resolution limit 
of the simulation. 

Unlike type II SNe and AGB stars, type la SNe are a result of 
binary stellar evolution, and consequently somewhat complicated 
to model. Currently there are two major theories for the progenitors 
of type la SNe. The most common is the single degenerate model. 
In this model, a white dwarf experiences enough mass transfer via 
a main sequence or giant companion to bring its mass above the 
Chandrasekhar limit, inducing an explosion. According to the dou- 
ble degenerate model, two white dwarfs merge after a long period 
of binary evolution. In order to predict the type la SN rate from 
a stellar population, one must consider a wide range of relatively 
uncertain processes (e.g. binary stellar evolution) and poorly con- 
strained parameters (e.g. binary fraction, binary separation, binary 
mass function). 

Although each SNIa is thought to inject a similar amount of 
kinetic energy as a type II SN, they are thought to be much less 
efficient at driving galactic outflows, mainly because the energy of 
a stellar population is released over billions of years rather than the 
tens of millions of years over which core collapse SNe release their 
energy. SNIa may, however, dominate the stellar energy feedback in 
galaxies with very low specific star formation rates. We therefore 
do include energy feedback by SNIa, which distribute in thermal 
form among the SPH neighbors of star particles. 

These two types of supernovae also have different chemical 
signatures. Type II SNe produce copious amounts of oxygen and 
so-called 'alpha elements' (neon, magnesium, silicon, etc.). These 
elements are primarily a result of a-capture reactions and, in addi- 
tion to oxyg en, make up the bulk of t he metal mass ejected from 
type II SNe dWooslev~ & Weaver 1995). What is known about type 
la SNe (SNIa) is that they are a major source of iron in the Universe 
and that a large fraction of them involve relatively old (> 10 9 yr) 
stars (see Greggio & Renzini 1983). Because of the time difference 
between the release of a elements by type II SNe and iron by SN 
la, the a/Fe ratio can provide information about the time since the 
last starburst, at least for the case of closed box models. 
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Figure Al. Stellar lifetime as a function of initial mass for various metal- 
licities (metal mass fractions). Solid lines show lifetimes as given by 
Portinari et al. 1 1998), while the dot-dashed line shows calculations from 
Padovani & MatteuccJ 1 19931) . who compiled results from the literature for 



a unspecified meta l licity, and the dashed line shows the lifetimes given by 
Maede r & Mevnej ll989h . who assumed solar abundances. Lifetime is a 
strongly decreasing function of mass, but it is only weakly dependent on 
metallicity. 



Al Stellar initial mass function 

The stellar IMF has long been a subject of debate. ISalpeteil il955h 
made the first attempt to derive an IMF from local stellar counts 
and his formulation is still widely used today. As star formation 
is still not understood well enough to predict the IMF from first 
principles, a wide ra nge of IMFs are available in the literature (e.g., 
Roma no et al.ll2005h . Although there is consensus about the fact 
that the IMF is less steep than Salpeter below 1 Mq , there are many 
uncertainties. In addition, the IMF may depend on redshift or on 
properties of the environment such as metallicity, gas pressure, or 
the local radiation field. 

For our purposes, the particularly interesting IMFs are the 
Salpeter IMF - because it serves as a refere nce mode l and b ecause 
it is still used in many simulations - and the Chabrier; (2003) IMF - 
because it is an example of an IMF with a low mass turnover that 
fits the observations much better than the Salpeter IMF. In this work 
we will assume a Chabrier IMF. 

The Salpeter mass function takes a very simple power-law 
form, 



$(M) 



dN 
dM 



AM' 



where the normalization constant A is chosen such that: 



M$(M)dM = 1 M e 



(Al) 



(A2) 



While the Chabrier IMF has the advantage that it does not overpro- 
duce low mass stars, it is slightly more complicated: 



M$(M) = 



Ae -dogM- 
BM~ 13 



log M y /la* 



if M < 1M 
if M > 1M 



(A3) 



where M c = 0.079, a — 0.69, and the coefficients A and B are set 
by requiring continuity at 1 Mq and by the normalization condition 
dA2t . Although the shape of the IMF above 1 Mq is very similar 
to Salpeter, the lognormal decrease at the low mass end results in a 
much lower stellar mass-to-light ratio. 

Besides the shape of the IMF, we also need to specify the mass 



limits. The lower limit is defined by the hydrogen burning limit of a 
star (typically cited as 0.08 M Q ; e.g. lKippenhahn & Weigertl 19941) . 
while the upper limit is significantly more uncertain, and the value 
most often found in the literature of 100 Mq roughly reflects the 
current upper limit of observed stars. For consistency with previous 
studies, we choose mass limits of 0.1 Mq and 100 Mq. 



A2 Stellar lifetimes 

Although a star can be active (via accretion) for an indefinite period 
of time, it is useful to define its stellar 'lifetime' as the time a star 
takes to move from the zero age main sequence up the giant branch 
and through any subsequent giant evolution. Using this definition, 
most stellar mass loss occurs at the end of the star's lifetime during 
the AGB/SN phase. 

There is no easy way to determine the lifetime function for a 
stellar population observationally. As a result, the published life- 
times _are_^uncdrjnal_ife to_tfie re sults of stellar evolution m odels 
(e.g.jMaeder & Mevnetll 19891 and lPadovani & Matteuccil 19931 : see 



iRomano et al. 1 20051 for an overview). 

Portin ari et al .1 J 19981) have published metallicity-dependent 
lifetimes from their stellar evolution calculations in tabular form 



and we show thei r results in figure 



used lifetimes of 



All together with the widely 



Maeder & Mevnet ( 1989), who assumed solar 



abundances, and Padovani & Matteuccj| 1 19931) . who presented a 



compilation taken from the literature but did not specify the as- 
sumed metallicity. The different lifetime sets agree well at high 
mass, but differ b y nearly an order of magnitud e for stars of a few 
solar masses with lPadovani & Matteuccil d 1993b giving systemati- 
cally lower values. 

Lifetime is a strongly decreasing function of mass and only 
weakly dependent on metallicity. Although metallicity seems to 
make little difference, it is attractive to be consistent in treating 
stellar evolution (in that we use yields based on the same calcula- 
tions). We thus emplo y the metallicity-dependent lifetime tables of 
|Portinariet"aH ( ll998l) . 



A3 Stellar yields 

The mass ejected in each element must be obtained from stellar 
evolution and (explosive) nucleosynthesis calculations. Differing 
treatments of physics can often lead to drastically different results 
in yield calculations. For instance, is it largely unknown what deter- 
mines the convection boundaries and hot-bottom burning in AGB 
stars, l eading to discrepanci es between yields setsF^l as we will see 
below. IHirschi et al.l {2005) found that adding rotation to type II 
SN models changes the yields by factors ranging from two to an 
order of magnitude. These problems are most evident in the fact 
that the yields of individual elements are often rescaled in order to 
get chemical evolution models to match observations. In fact, some 
authors (e.g. iFrancois et al. 2004) use observations and chemical 
evolution models to derive stellar yields for a number of elements, 
and the results show marked differences from yields derived from 
nucleosynthetic calculations. 

Another important point to make sure of when using yields 
from various processes, is that the theoretical assumptions made in 
the models match. For instance, if one uses AGB yields for masses 
up to 8 Mq, then one must take care to use type II SNe yields that 



10 In this section we use the term yield in the generic sense, referring to 
stellar ejecta. We will define it more rigorously in appendix !4. II 
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Table Al. AGB yield references and grids 



Reference 


Initial stellar mass (Mq ) 




Metallicity (metal mass fraction) 


van den Hoek & Groenewegen (1997) 


[0.8, 0.9, 1, 1.3, 1.5, 1.7, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 


8] 


[0.001, 0.004, 0.008, 0.02, 0.04] 


Marigo(2001) 


(various) [0.85 - 5] 




[0.004, 0.008, 0.019] 


Izzard et al. (2004) 


[0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 


8] 


[0.0001,0.004, 0.008,0.02] 




10" J io~ 

Stellar Initial Metal Mass Fraction 




-0.5 



10" J 10" 
Stellar Initial Metal Mass Fraction 



Figu 

ot van den Hoe k & Groenewegen] Jl997l) (black), iMarigcl feOOll) (red), or Izzard et al. (2004) (blue). Shown are the abundances of helium (left, solid), carbon 
(left, dashed), nitrogen (right, solid) and oxygen (right, dashed) relative to hydrogen in solar units. The calculations assume a Chabrier IMF and integrate the 
yields over the stellar initial mass range [0.8, 6] Mq. The different yield sets agree around solar metallicity, but differ significantly for lower metallicities. 



begin higher than 8 Mq dMarigol200ll) . The problem is more often 
than not in reverse, however, since most yield sets for massive stars 
are tabulated down to ~ 10 Mq, but most intermediate mass yield 
sets only go up to at most 8 Mq, leaving it up to the user to decide 
what to do for the transition masses. Ideally, one would like to use a 
consistent set of yields, in the sense that the same stellar evolution 
model is used for both intermediate mass stars and the progenitors 
of core collapse SNe. 



A3.1 AGB stars 

The asymptotic giant branch phase of stellar evolution occurs in in- 
termediate mass stars near the very end of their lifetimes. During 
this phase the envelope of the star puffs up and is eventually shed, 
causing the star to lose up to 60 % of its mass. Prior to the AGB 
phase, material in the core (where most of the heavy elements re- 
side) is dredged up into the envelope via convection. As a result, 

the ejecta are particularly rich in carbo n and nitrogen. 

Building on pion eering work by llben & Truranl and 
iRenzini & Voli 1 198 lb. various groups have published AGB yields 



(e.g. Forestini & Charbonnel 1997; van den Hoek & Groenewegen 
19971; IMarigol 1200 ll ; IChieffi et al.l 1200 ll ; iKarakas et alj |2002| ; 



Izzard et al 



2004T) . Table |A]J outlines the extent of the 
res olution (in mass and meta ll icity) o f the AGB yields 
of Ivan den Hoek & Groenewegen! J 19971) . IMarigol J200ll) and 
Izzard et al. I ( 120041) . which are some of the most complete sets for 
our purposes. These yields are compared in Fig. IA2I which shows 
the abundance relative to hydrogen, in solar unita 11 L of various el- 



We use the notation [X/H] = log(X/H) - log(X/H) Q . 



ements in the ejecta as a function of stellar metallicity. These cal- 
culations are for an SSP with a Chabrier IMF and the mass range 
[0.8, 6] Mq at time t = oo. The yields agree very well at solar 
metallicity, and for the case of helium, this agreement extends to 
lower metallicities. However, for nitrogen, oxygen, and particularly 
carbon, different yields sets give very different results at low metal- 
licities. 

We show in Fig. IA3I for each yield set, the integrated frac- 
tion of the initial SSP mass ejected by stars in the mass range 
[0.8, 6] Mq at time t = oo, normalized to the total initial stellar 
mass over the range [0.1, 100] Mq. The ejected mass fractions are 
very similar for the different yield sets. 

In this work we use the yields of Marigo] d200ll) . Although 
these only go up to 5 Mq, they form a complete set with the SN 
Type II yields of IPortinari et al.1 d 19981) since they are both based 
on the Padova evolutionary tracks. Indeed, there are very few yield 
pairings that form a consistent set across the full mass range. 



A3. 2 Type II Supernovae 

Although major advances have been made in modelling type 
II SNe, theorists have yet to come up with a model that self- 
consistently explodes. As a result, the yields are very sensitive to 
the location of the mass cut, i.e., the explosion radius. The material 
interior to this radius is assumed to end up in the stellar remnant, 
while the mass exterior to this radius is assumed to be ejected in 
the explosion. Unfortunately, the mass cut is very uncertain, thus 
creating another degree of freedom. Although uncertainties in the 
explosion radius do not affect the oxygen yields much, the yields 
of heavier elements such as iron are very sensitive to this choice. 
A large number of groups have published tables of type II SN 
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Stellar Initial Metal Mass Fraction 

Figure A3. Total fraction of the initial mass of an SSP ejected by AGB stars 
at time t = oo as a function of initial stellar metal mass fract ion, assuming 
the yields of van den Hoek & G roenewegen J 19971) (fcfacfcl. lMarigol fcOOll) 
(red), or Izzard et al. 1 2004) (blue). The calculations assume a Chabrier IMF 
and integrate the yields over a stellar initial mass range [0.8, 6] Mq, but 
the fractions are normalized to the mass range [0.1, 100] Mq. Different 
yield sets predict similar ejected mass fractions. The ejected mass fraction 
is insensitive to metallicity. 



yields (e.g..lMaedeill992l:IWooslev & Weaver! 1 9951: IPortinari et al 
19981: iRauscher et alj 120021 : iHirschi et al.l 120051: iNomoto et al 
2006). Although the more recent models tend to be more sophisti- 



cated, for example including the effec ts of improved stel lar and nu- 
clear physics and rotation, the Wooslev & Weaver ( 1995) yields are 
still the most widely used. IPortinari et alj |T998) combine models 
of pre-SN s tellar mass loss with the nu cleosynthesis explosion cal- 
culations o f IWooslev & Weaver! d 19951) . who ignored pre-SN mass 
loss. To do this they link the carbon-oxygen core mass p redicted by 
their models with those of IWooslev & Weavetj d 19951) . They find 
low metallicity, massive stars to have very inefficient mass loss for 
these large core masses (for which there are no type II SN yields), 
they only consider mass loss by winds and assume the rest of the 
star collapses directly into a black hole. 

Tabl e IA2 | and Figure IA4I c ompare the SN typ e II 
yields of IWooslev & Weaveil d 19951) . IPortinari et ail dl998l) . and 
Chieffi & Limongi (2Q0J). Fi gure IA4I shows the abundance rela- 
tive to hydrogen, in solar units, of various elements in the ejecta as 
a function of stellar metallicity. The elements shown are those that, 
together with hyd rogen, dominate the radiative cooling of (photo- 
)ionized plasmas dWiersma et al .1 12009). These calculations are for 
an SSP with a Chabrier IMF in the mass range [8, 40] M Q at time 
t = 00. It is clear that the yields of elements that are produced by 
type II SNe depend only weakly on metallicity, unless the metal- 
licity is supersolar. However, except for very low metallicities, the 
nitrogen abundance in the ejecta is proportional to the stellar metal- 
licity, indicating that it is simply passing through rather than being 
produced. 

The different yield sets agree well for helium, carbon, and ni- 
trogen, but there a re large differences fo r he avier elements. The 
differe nce between IPortinari et al. I dl998h and IWooslev & Weaveil 
d 19951) is due to the fact that the former add their stellar evo- 
lution to the IWooslev & Weaveil d 19951) nucleosynt h esis c alcu- 
lations. The difference between Chieffi & Limongil d2004l) and 
Woosl ev" & Weaveil d 19951) is caused mostly by the difference in 
the assumed mass cut. 

The yields presented in IWooslev & WeaveJ d 19951) consider 



the state of the ejecta 10 5 s after the explosion. Because a number of 
isotopes have rather short decay times, it is customary to consider 
- as most recent yield sets have do ne - the state of the ejec ta at a 
much later time (10 8 s in the case o f lChieffi & Limongil20 04). This 
is esp ecially important for 5 6 Ni, which decay s rapidly into 56 Fe. 
When IPortinari et alj d 19981) incorporate the Woos lev & Weaveil 
d 19951) yields into their calculations, they simply add the 56 Ni to 
the 56 Fe. We have added an additional magenta l ine to figure lA4l 
to show the iron yield from I Wooslev & Weaved dl995h with this 
adjustment. This agrees much better with the other two yield sets. 

Figure fA5] shows the integrated fraction of the initial mass of 
an SSP that is ejected by type II SNe at time t = 00 as a function 
of metallicity. These calculations again assume a Chabrier IMF and 
integrate the yields over the mass range [8, 40] Mq, but normalize 
the mass fraction to the mass range [0.1, 100] Mq. The ejected 
mass is insensitive to metallicity and the three yield sets are in ex- 
cellent agreement. 

We use the yields of IPortinari et alj £l998) since they in- 
clude mass loss from massive stars and because they form a self- 
consistent set toge ther with the ste llar lifetimes of these authors and 
the AGB yields o f lMarigol d200lh . Their models allow for the pos- 
sibility of electron capture supernova for intermediate mass stars 
(7 — 9 Mq). Unfortunately, nucleosynthesis calculations for thes e 
star s have only recently b een performed (e.g. JWanaio et al]|2008l) . 
The lPortinari et al] d 19981) yields therefore only consider the shed- 
ding of the envelop e and t he m ass loss up to that stage. Wh ile both 
IChieffi & Limongil d2004l) and IWooslev & Weaveil d 1995b evolve 
their models until the point of supernova and then begin super- 
nova calculations, they do not include mass loss in their yield tables 
(since the goal of their studies was to investigate explosive nucle- 
osynthesis). 

Following the recommendation of L. Portinari (private com- 
munication), we have adjusted the massive star yields by the fac- 
tors inferred from comparisons of che modynamical models and 
current Galactic abundances. Following IPortinari et aD{l998), wc 
have multiplied the type II SN yields of C, Mg, and Fe for all 
masses and metallicities by factors of 0.5, 2, and 0.5, respec- 
tively. Note that these factors were not included in Fig. IA4I The 
adjustments to Mg and Fe can be justified due to uncertainties 
in the explosive nucleosynthesis models by IWooslev & Weaver 
d 19951) as discussed by a number of authors (e.g. iTimmes et al.l 
1995l:lLindner et al.ll999l:ICarigil2000l:lGoswami & Prantzosl200l 
Carigil 12003k iGavilan et alj 120051 : iNvkvtvuk & Misheninal I200T 



These a djustments have also been incorpo r ated into other stud- 
ies (e.g. ilianget all 1200 ll ; Ilia et alj [20021 ; IPortinari et al] |2004 
iNagashima et alj|2005h . These ad hoc adjustments reveal just how 
uncertain the yields are. 



A3. 3 Type la Supernovae 

Despite the fact that the progenitors of type la SNe are still un- 
certain, yield calculations can be made based on an explosion of 
a Chandrasekhar mass carbon-oxygen dwarf. While most models 
assume the progenitor to be a binary system, the yields usually 
consider just the explosion of the white dwarf itself, ignoring the 
companion completely. 

It is instructive to anticipate a few features of these calcula- 
tions. First, since a type la SN explosion occurs as soon as the 
compact object reaches the Chandrasekhar mass, the yields from 
this process will be independent of the initial stellar mass. Sec- 
ond, by the time a star gets to the white dwarf phase, it is almost 
completely composed of carbon, nitrogen and oxygen, and while 
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Table A2. SN type II yield references and grids 



Reference 


Initial stellar mass (Mq) 


Metallicity (metal mass fraction) 


Woosley & Weaver (1995) 
Portinariet al. (1998) 
Chieffi & Limongi (2004) 


[11, 12, 13, 15, 18, 19,20, 22, 25, 30, 35,40] 
[6, 7, 9, 12, 15, 20, 30, 60, 100, 120, 150, 200, 300, 500, 1000] 
[13, 15, 20, 25, 30,35] 


[0.0, 0.000002, 0.0002, 0.002, 0.02] 
[0.0004, 0.004, 0.008, 0.02, 0.05] 
[0.0, 0.000001, 0.0001, 0.001, 0.006, 0.02] 




Figur e A4. Composition of t he integr ated SN type II ejecta of an SS P at time t = oo as a funct ion of its initial stellar metal mass fraction, assuming the yields 
of lWooslev & WeaveJ 1 19951) (fefad:). |Portinari et aljfl998l) (red), or lChieffi & Limongi 120041) (blue). Shown are the abundances of helium, carbon, nitrogen, 
oxygen, neon, magnesium, silicon, sulphur, calcium, and iron, as indicated in the legends. The calculations assume a Chabrier IMF and integrate the yields 
o ver the stellar initial mass ra nge [8, 40] Mq. The (magenta) dashed line in the lower right panel indicates the result of transferring the 5B Ni to the Iron yield 
of Wooslev & We aver! 1 19951) . While the different yield sets agree well for the lightest elements (upper panels), there are large differences for elements heavier 
than nitrogen. 



the balance between these elements may depend on initial mass or 
composition, the models assume such a dependence is insignificant. 



Several research groups have published yields from type la 
SNe. Their calculations range fr om one- to three-dimensional 
calculations I Travagl io et al . 2004). The standard to which most 
results are com pared is the spher i cally s ymmetric calcu lation 
dubbed "W7" jNomoto et al J 1 1984 1 19971: llwamoto et all [l999l ; 
iBrachwitz et al.lr2000l : iThielemann et al.1120031) . We use the latest 
inc arnation (at the time of c ode development) of the W7 model, 
i.e. . IThielemann et alj d2003l) . 



A4 SN type la rates 

The recipes for the release of mass by AGB stars and SNII are quite 
simple because these processes are direct results of stars reaching 
the ends of their lifetimes. Hence, one merely needs to combine the 
IMF with the stellar lifetime function and the yields and integrate 
over the time interval specified. The recipe for SNIa is, however, 
somewhat more complex. 

Two SNIa channels a re thought to be most plausible (see, e.g., 
IPodsiadlowski et al]|2008l for a review). According to the single 
degenerate scenario, SNIa result when accretion onto a white dwarf 
by a binary companion pushes the former over the Chandrasekhar 
mass. The double degenerate scenario, on the other hand, involves 
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1CT = 10 1CT J 1CT Z 10"' 

Stellar Initial Metal Mass Fraction 

Figure A5. Total fraction of the initial mass of an SSP ejected by SNe 
type II at time t = oo as a function of init i al stel lar metal mass fraction, 
assum ing the yie lds of IWooslev & Weaveil jl995l) (black), [Portinari et al. 
1 19981) (red), or IChieffi & Limongil 120041) (blue). The calculations as- 
sume a Chabrier IMF and integrate the yields over the stellar initial mass 
range [8,40] Mq, but the fractions are normalized tot he mass range 
[0.1, 100] Mq. The ejected mass fraction is insensitive to metallicity and 
different yields sets predict very similar results. 



the merger of two white dwarfs, thus putting the merger product 
over the Chandrasekhar mass. 

The SNIa rate of an SSP as a function of its age can be deter- 
mined using late age stellar and binary e volution theory. The formu- 
lation that is used most often is given bv lGreggio & Renzinil J 1983h . 
This formulation requires knowledge of the distribution of sec- 
ondary masses in binaries and makes some assumptions about bi- 
nary evolution. These are poorly constrained and may therefore be 
subject to large uncertainties. Hence, it is attractive to take an alter- 
native, more empirical approach. Guided by observations of SNIa 
rates, we can simplify the prescription by specifying a functional 
form for the empirical delay time function £(t), which gives the 
SNIa rate as a function of the age of an SSP, normalized such that 
J"°° £(t)dt = 1. The parameters of this function can then be deter- 
mined by fitting to observ ations of SNIa rates (e.g.. lBarris & Tonrvl 
l2006l : lForster et al.ll2006h . 

The number of SNIa explosions per unit stellar mass in a time 
step At for a given SSP is then 



N sm *(t;t + At) = v 



t+At 



itp')dt', 



(A4) 



where v is the number of SNe per unit formed stellar mass that 
will ever occur. While this approach is attractive since it does not 
require us to make any assumptions regarding the progenitors, we 
will employ yields that correspond to a scenario involving at least 
one white dwarf so the SNIa rate should take that into account 
(i.e., we should not have any SNIa before the w hite dwarfs have 
evolve d). We therefore take an approach similar to lMannucci et al.l 
(2006) and write: 



N S m*(t;t + At) 



t+At 



Mt')Z(t')dt', 



(A5) 



where a is a normalization parameter and / w< j (t) is the number of 
stars that have evolved into white dwarfs up until time t (i.e., the 
age of the SSP) per unit initial stellar mass: 



/wd(t) = 



I" m wdhigh 



if M z (t) > m w dhigh 



$(M)dM otherwise 



(A6) 



' m SNIalow(*) 



Here m w dhigh and m w diow are the maximum and minimum white 
dwarf masses respectively, Mz (t) is the inversj^l of the lifetime 
function tz(M), and 



m-SNiaiow(i) = max(M z (t),m wd i ow ) 



(A7) 



Note that the shape of the SNIa rate as a function of time differs 
between equations dA4b and JA5b . 

It is thought that stars with main sequence masses between 
3 Mq and 8 Mq evolve into a SNIa progenitor white dwarf. Note 
that while our type II SN yields range to m asses as low as 6 Mq , 
the models used in the yield calculations of IPortinari et al.l i 19981) 
for stars of 6 and 7 Mq do not incorporate any explosive nucle- 
osynthesis. They note that stars of these masses could explode as 
electron capture SNe or evolve to a thermally-pulsing AGB phase, 
thus simply shedding their envelopes. 

The form of the delay function £ h as generated particular in- 
terest recently (e.g.. iDahlen et al.ll2004f) . The two types of delay 
functions that are most often considered, the so-called e-folding 
delay and Gaussian delay functions, are shown in figure lAol The 
e-folding delay function takes the form: 



Tla 



(A8) 



where ri a is the characteristic delay time. This delay approximates 
predictions made for the single degenerate scenario via population 
synthesis models. 

The Gaussi an delay function w as motivated by high redshift 
observations by Dahl en et alJd2004T) . which show a marked decline 
in the SNIa rate beyond z = 1. It takes the form: 



V27I 



(A9) 



where a — 0.2ri a for 'narrow' distributions and a — 0.5ri a 
for 'wide' distributions. Note that the integral of this particular 
function (J f(t)dt is only normalised to unity for a <C Ti a . 
Gaussian delay functions often feature long characteristic times 
(ria = 4 Gyr) in order to fit the data. Strengths and weaknesses 
exist for both delay functions, and are discussed in the references 
cited above. Unfortunately, the choice of distribution (including all 
the ensuing parameters) is rather poorly constrained. In particular, 
both the cosmic star formation rate and the type la rate contain large 
uncertainties beyond a redshift of 1. For in stance, the type la rate 
quoted at z = 1.6 bv lDahlenetai] ( l2004l) is based on only two 

detection events. 

In addition to the cosmic SNIa rate, iMannucci et al. I d2006l) 
attempt to use other observations to constrain the shape of the de- 
lay function. They cite the dependenc e of the SNIa rate on galaxy 
colour observed in the local universe (Mannucci et al J2005I) as well 
as a d ependence on radio loudness observed by iDella Valle et all 
d2005h . Their analysis suggests that the delay time function could 
hav e two components, a "prompt" mode and a "tardy" mode (see 
also lScannap ieco & Bildsten 2005). These modes may have phys- 
ical counterparts with the two progenitor channels, but this is still 
uncertain. Using such a delay function could improve our fits to 
the observations marginally, but would introduce another unknown 

12 The lifetime function is invertible because it is a monotonic function of 
mass for a fixed metallicity. 
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Figure A6. Volumetric Type la supernova rate as a function of redshift 
from our simulations (black) using either an e-folding {solid) or a Gaus- 
sian (dashed) delay function. Also shown are approximate fits using a 
"standard formulation" of the cosmic SNIa rate, calculated using equa- 
tion )A4t and the star formation history predicted by the simulations and 
a Chabrier IMF (red ) . The data points corr espond to observa t ions r eported 
by ICappellaro etaT] Jl999l) (filled c ircle), iMadgwick et alj 120031) (open 
square ), iBland J2004I) (fi lled square). lHardirj feOOoF (open diamond), Neill 
feOOd) (filled diamond). iTonrvl 120031) (open black triangle) | PairJ j2002l) 
(filled t riangleVlPahlen et alj J2004I) fcrosses). lBarris & Tonrvl KOOd) (open 
circles ), |PoznanskTet*aTn2007l) (upright triangles), and iKuznetsova et alj 
1 2008) (left pointing triangles). 



parameter (the ratio of contributions of the modes). Moreover, the 
e-folding delay function also includes significant prompt and late 
contributions. We therefore chose to use an e-folding delay func- 
tion with ria = 2 Gyr, but have also performed one simulation with 
a Gaussian delay function using n a = 3.3 Gyr and a = 0.66 Gyr. 
These parameter values were chosen to roughly agree with the con- 
straints mentioned above. The coefficient a appearing in equation 
l |A5t was chosen to roughly match observations of the cosmic SNIa 
rate. 

We plot the current measurements of the cosmic SNIa rate 
in figure |A6l We have self-consistently adjusted all data points to 
our cosmology of choictEl {h,D. m ,D. A ) = (0.73,0.238,0.762). 
Note that the different SNIa measurements are strongly discrepant, 
indicating that the statistical errors have been underestimated 
and/or that some rates suffer from systematic errors. 

The solid, black curve in Figure lA6l shows the evolution of the 
SNIa rate in our L100N512 simulation, which used the e-folding 
time delay function and a = 0.01. Also shown (dashed, black 
curve) is the predicted type la rate for another simulation that is 
identical to L100N512 except for the fact that it uses the Gaus- 
sian delay function and a — 0.0069. The e-folding and Gaussian 
models result in reduced x 2 values of 2.6 and 2.2, respectively. We 
consider this acceptable since the data are internally inconsistent. 
Note that the predicted SNIa rates agree better with the more recent 
measurements. 



13 Correcting for cosmology is important since observations are taken 
over a volume in co-moving space. For instance, the highest redshif t 
point of lPoznanski et alj 120071) is greater than that of lDahlen et alj J2004I) 
when using their c osmology, but when converting to our cosmology the 
iDahlen et al]<2004l) point is reduced by a greater amo unt since it was taken 
over a larger redshift rang e. Thus, our plot shows the lDahlen et alj J2004I) 
point to be greater than the Poznanski et al. 1 2007) point. 



While preparing this publication, we encountered a bug in the 
code that resulted in m w diow being set to zero in equation {AT}. The 
net result is a shift in the effective delay time to lower redshifts and 
an increase in the effective value of a. Because this error compli- 
cates the interpretation of our parameter values and because most 
of the literature uses equation dA4t rather than ( |A5I > we have tried 
to approximately match the SNIa rates predicted by our L100N512 
simulations using equation JA4t . taking the star formation history 
predicted by the simulations as input and using the same (Chabrier) 
IMF as was used in the simulation. We expect the effects of this bug 
to be nominal since our rates still pass comfortably though the ob- 
servations. Moreover, a simulation with a Gaussian delay function 
(which we will present in a later publication) showed negligible 
difference to the e-folding model in nearly all respects (e.g., star 
formation history, distribution of metals). The difference between 
the type la rates is much greater than the difference between the 
bugged and unbugged versions. 

The matching "standard formulation" SNIa rates are shown as 
the red curves in figure IA61 The e-folding delay function (solid, red 
curve) uses rt a = 3 Gyr while the Gaussian delay function (dashed, 
red curve) uses n a = 3.3 Gyr, a = 0.66 Gyr (which is identical to 
the original delay function used in combination with equation I A5b . 
When using equation dA4t it is useful to parametrize the normal- 
ization in terms of r\, the fraction of white dwarfs that eventually 
explode as SNIa, or the Type la efficiency, which is related to v, the 
number of SNIa per unit formed stellar mass via 

r\ I $(M)dM = v I A/$(M)dM. (A10) 

J3 JO.l 

The red curves in figure |A6l correspond to 77 = 2.50 % for the e- 
folding model and r\ — 2.56 % for the Gaussian model. Note that 
these efficiencies correspond to our Chabrier IMF, whereas much 
of the literature quotes efficiencies for a Salpeter IMF. 

We end the discussion by noting that we match another con- 
straint in that we find iron abundances in our galaxies to be roughly 
solar at z — 0, as we will show elsewhere. 



APPENDIX B: VARYING THE SIZE OF THE 
SIMULATION BOX 

In this appendix we use the suite of cosmological simulations listed 
in Table [2] to test for convergence with the size of our simulation 
box. The size of the simulation box is important because it deter- 
mines what kind of objects can form. Rare, large-scale structures 
(both high- and low-density peaks) can obviously only be sampled 
correctly if they are much smaller than the size of the box. More- 
over, since Fourier components of the density field only evolve in- 
dependently in the linear regime, the simulation box must be large 
compared with the scale that has last gone non-linear. Because this 
scale increases with time, larger boxes are required at lower red- 
shifts. To isolate the effect of the size of the simulation volume, we 
will compare simulations that use different box sizes, while holding 
the resolution constant. 

For this test we use two sets of simulations that dif- 
fer in terms of their resolution. The first set consists of 
L100N512, L050N256, and L025N128 and the second set com- 
prises L025N512, L012N256, and L006N128. Within each set the 
size of the box is thus decreased by factors of two and four, respec- 
tively. The second set uses a particle mass that is 64 times smaller 
than the first set, but these higher resolution runs were not contin- 
ued down to z = 0. 
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Figure Bl. Dependence on the simulation box size of the evolution of the fractions of the metal mass in various components. The curves in the top panels 
correspond to the L100N512 (black), L050N256 (red), and the L025N128 (blue) simulations. The curves in the bottom panels are for the L025N512 (black), 
L012N256 (red), and the L006N128 (blue) simulations. The line styles are identical to those shown in Fig. [7] The parts of the curves corresponding to redshifts 
for which the total metal mass is smaller than 10 — 6 of the total baryonic mass have been omitted because they become very noisy. Except for the ICM 
(T > 10 7 K), the 50 h' 1 Mpc and 12.5 h' 1 Mpc boxes have nearly converged for z < 2 and z > 2, respectively. 




Figure B2. Dependence on the simulation box size of the evolution of the metallicities of various components. The colors and line styles are identical do those 
used in Fig. lBll The 25 h~ 1 Mpc and 12.5 h" 1 Mpc boxes have converged for z < 2 and z > 2, respectively. 



The top and bottom rows of fi sure IbTI show the evolution of 
the metal mass fractions in the various components for the low and 
high-resolution sets, respectively. The panels are analogous to those 
shown in Fig. [7] In particular, the black curves in the top panels, 
which correspond to L100N512, are identical to those shown in 
Fig-El The top panels demonstrate that, except for the T > 10 5 K 
gas (right panel), the results have already nearly converged for 
the 25 h~ x Mpc box. For the WHIM (right; solid) we require a 
50 h^ 1 Mpc box, although even the 25 h^ 1 Mpc box seems to be 
nearly converged for z < 2. For the ICM (right; dashed), how- 
ever, even the 50 h^ 1 Mpc box is too small. The bottom panels 
demonstrate that while a 12.5 h~ x Mpc box is sufficient for z > 2 
(except for the ICM, which accounts for too few metals to be vis- 



ible in the plot, and the photo-ionized IGM), boxes as small as 
6.125 h~ l Mpc give spurious results for z < 6. The fact that we 
require larger boxes for hot gas is not surprising because objects 
with high virial temperatures are rare and large. 

The evolution of the metallicities of the different components 
in the two sets of simulations is shown in Fig. IB2I Interestingly, 
except for the ICM, the results appear to have converged already 
for our smallest simulation boxes (although this is not quite true 
for the stellar metallicity in L006N128). Together with Fig. IBTI this 
suggests that it is easier to obtain convergence for the metallicity 
than for the metal mass fraction. This in turn implies that the lack of 
convergence of the metal mass fractions for the smallest boxes was 
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Figure B3. Dependence on the size of the simulation box of the probability density function, weighted by metal mass, of the gas density (left), temperature 
(centre-left), gas metallicity (centre-right), and stellar metallicity (right) at z = (top) and z = 2 (bottom). The colors and line styles are identical to 
those used in Fig. IB II The vertical, dotted lines in the left panels indicate the threshold for star formation. Note that star-forming gas was excluded from the 
temperature panels. Except for the extremes of the distributions, the 50 h^ 1 Mpc and 12.5 h -1 Mpc boxes have nearly converged for z < 2 and z > 2, 
respectively. 



caused by non-convergence of the baryonic mass fraction rather 
than the metallicity. 

However, part of the difference in appearance can be ac- 
counted for by the fact that the y-axis spans six decades in Fig. lB2l 
but only two decades in Fig. IB II Closer inspection reveals that 
while this explanation may be viable for the stars, it does not hold 
for the gaseous components. Furthermore, the axis ranges differ for 
a reason: they show the range of interest. While we generally do not 
care whether the metal mass fraction in a particular component is 
1CF 3 or 1CP 2 when drawing up a census of metals (both fractions 
are negligible), it is very interesting to know whether the diffuse 
IGM has a metallicity of 10 -3 or 10~ 2 solar because the differ- 
ence is measurable and because it has important consequences for 
enrichment scenarios. 

Finally, we plot the probability density functions of the gas 
density, temperature and metallicity as well as the stellar metallic- 
ity, all weighted by metal mass, in Fig. IB3I The results confirm the 
conclusions drawn from the previous two figures. This figure illus- 
trates clearly that differences first show up at the extremes of the 
distribution, which is not surprising since those correspond to rare 
objects. 

Summarizing, a 50 Mpc box is sufficient to obtain a con- 
verged picture of the cosmic metal mass fractions and metallicities 
of most components down to z — 0. Somewhat smaller boxes may 
suffice if one is only interested in higher redshifts (except for the 
ICM, 12.5 bT x Mpc is sufficient for z > 2) or in low-mass ob- 
jects. We caution the reader that it is possible that larger box sizes 
are needed if other aspects of the cosmic metal distribution are in- 
vestigated (e.g., clustering strengths of metal absorption lines). 



APPENDIX C: VARYING THE RESOLUTION 

While the size of the simulation box mainly determines the types 
of structures that can form, numerical resolution can even have a 
large effect on common objects. For example, simulations that do 
not resolve the Jeans scales may underestimate the fraction of the 
mass in collapsed structures and hence the total amount of metals 
produced. Moreover, processes like metal mixing (see $5.2\ may 
depend on resolution. To isolate the effect of numerical resolution, 
it is important to hold the simulation box size constant (or to use a 
box sufficiently large for the result to have converged with respect 
to the size of the simulation box). 

Before showing the results of the convergence tests, it is use- 
ful to consider what to expect. The temperature of substantially 
overdensj* 4 ! gas does not drop much below 10 4 K in our simu- 
lations (see Fig. |9). In reality, gas at interstellar densities (nu > 
1CP 1 cm" 3 ) is sufficiently dense and self-shielded to form a cold 
(T <C 10 4 K), interstellar phase, allowing it to form stars dSchavd 
2004). However, our simulations impose an effective equation of 
state for gas with densities exceeding our star formation threshold 
of riH = 10 _1 cm -3 . For our equation of state, P oc p 4 ^ 3 , the 
Jeans mass is independent of the density. Hence, provided we re- 
solve the Jeans mass at the star formation threshold, we resolve it 
everywhere. The Jeans mass is given by 

«, „ WO'.-'Me/r (j^f" CO 

where / g is the local gas fraction. Hence, we do not expect con- 
vergence unless the gas particle mass m g <C 10 7 Mq. To achieve 

14 Gas with very low overdensities can have temperatures substantially be- 
low 10 4 K due to adiabatic expansion, but the Jeans scales corresponding 
to these low densities are nevertheless large. 
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convergence, a simulation will, however, also need to resolve the 
Jeans length Lj, which implies that the maximum, proper gravi- 
tational softening, e P ro P , must be small compared with the Jeans 
Length 

Note that since Lj scales as Lj oc p -1 ' 3 for our equation of state, 
the softening scale will always exceed Lj for sufficiently large den- 
sities. However, since the Jeans mass does not decrease with den- 
sity, runaway collapse is not expected for star-forming gas. Com- 
paring the above equations with the particle mass and softening 
scales listed in Table [2] we see that while our highest resolution 
simulation (L025N512) marginally resolves the Jeans scales for 
/ g w 1, this is not the case for the simulations that go down to 
2 = 0, although L050N512 has m g ~ Mj and e pr0 p ~ Lj and is 
therefore not far off. 

However, resolving the Jeans scales post-reionization may not 
even be sufficient, given that the simulations start at 2 > z rc ion. 
Prior to reionization the gas temperature can be much lower in 
which case we have no hope of resolving the Jeans scales at 
the threshold for star formation. On the other hand, a photo- 
dissociating UV background may well prevent the collapse of 
haloes with virial temperatures smaller than 10 4 K. In any case, we 
expect the duration of this period to be comparatively short. Our 
simulations assume « fe ion = 9 and include a photo-dissociating 
background for z > 9 (see Sj2]l. Hence, in our simulations the for- 
mation of haloes with T v i r <§C 10 4 K will be efficiently suppressed. 

We test for convergence with resolution using two sets of sim- 
ulations, which use different box sizes. The first set consists of 
L100N256, L100N512, and L050N512 and the second set com- 
prises L025N128, L025N256, and L025N512. Within each set the 
mass (spatial) resolution is thus increased by factors of 8 and 64 
(2 and 4), respectively. Note that L050N512, the highest resolution 
run of the first set, uses a 50 /i -1 Mpc box whereas the other runs 
of this set used a 100 h~ x Mpc box. However, as we have shown 
in appendix [B] 50 ft -1 Mpc is sufficiently large to give converged 
results for all phases but the ICM. We do not show L100N128 be- 
cause most of its predictions are completely unreliable due to its 
extremely low resolution. The second set of simulations uses a box 
size of only 25 h~ x Mpc, but these higher resolution runs were not 
continued down to z = 0. 

Fi gures [CTl - IC 3 1 are analogous to fi gures IbTI — IB3I but show 
the effect of varying the resolution rather than the size of the sim- 
ulation box. The top row of Fig. IC 1 1 shows no evidence for full 
convergence, but the difference between L100N512 (black) and 
L050N512 (dark green) is small for z < 1. From the bottom row 
of panels we can see that L025N512, while certainly not fully con- 
verged, is similar to L025N256 for z < 4. As the resolution is 
increased, the fractions of the metal mass in stars and the WHIM 
tend to increase, while the fractions in non-star-forming gas with 
T < 10 5 K tend to decrease. 

It is not surprising that convergence is better at lower redshift, 
because the typical mass of star-forming haloes increases with de- 
creasing redshift which means it will be sampled with more parti- 
cles per halo. In particular, lower resolution simulations will start 
forming stars later and it takes some time for the resultant underes- 
timate of the metal mass to become negligible. 

Figure IC2l shows that it is much easier to get convergence for 
the metallicity of the different components than it is to get con- 
verged predictions for the metal mass fractions. As for convergence 
with box size, this partly reflects the fact that we do not require as 



much precision for the metallicities since the dynamic range of in- 
terest is much larger than it is for the metal mass fractions. The 
metallicities are nearly converged in the L050N512 run for z < 2 
and in the L025N512 for z < 4. Finally, Fig. lC3l confirms the find- 
ing from Fig. IC 1 I that increasing the resolution mainly moves met- 
als from cold-warm halo and ISM gas to the WHIM. In addition, 
Fig. IC3I shows that the high-density cut-off to the gas distribution 
is mainly set by resolution, which is expected because many parti- 
cles and small softening lengths are needed to sample the tail of the 
density distribution. The main conclusion we draw from this figure 
is, however, that the convergence is sufficiently good to draw inter- 
esting conclusions, particularly for the 25 /i -1 Mpc box at z = 2 
and the 50 h' 1 Mpc box at z = 0. 

Thus, the results of the convergence tests mostly agree with 
the expectations based on comparisons of the numerical resolu- 
tion with the Jeans scales. The resolution of the L100N512 run 
(ra 9 « 9 x 10 r ft" 1 Mq) is sufficient to draw qualitative conclu- 
sions regarding the mass fractions in various phases and for z < 2 
the L050N512 (m 9 » 1 x 10 7 h' 1 M Q ) may be close to conver- 
gence. Run L100N512 already has sufficient resolution to obtain 
interesting, quantitative predictions for the metallicities of the dif- 
ferent phases at z < 2. For z > 2 we need higher resolution. 
While the mass fractions are still not fully converged for L025N5 1 2 
(nig ~ 1 x 10 6 hT 1 Mg), its predictions for the metallicities are 
robust for z < 4. 
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Figure CI. Dependence on the numerical resolution of the evolution of the fractions of the metal mass in various components. The curves in the top panels 
correspond to the L050N512 (dark green), L100N512 (black), and the L050N256 (red) simulations. The curves in the bottom panels are for the L025N512 
(black), L025N256 (red), and the L025N128 (blue) simulations. The line styles are identical to those shown in Fig. [7] The parts of the curves corresponding 
to redshifts for which the total metal is smaller than 10~ 6 of the total baryonic mass have been omitted because they become very noisy. Convergence is 
relatively poor for z > 2 in the top row and for z > 4 in the bottom row. 




Figure C2. Dependence on numerical resolution of the evolution of the metallicities of various components. The colors and line styles are identical to those 
used in Fig. lCll Convergence is good for z < 2 in the top row and for 2 < z < 4 in the bottom row. 
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Figure C3. Dependence on numerical resolution of the probability density function, weighted by metal mass, of the gas density (left), temperature (centre-left), 
gas metallicity (centre-right), and stellar metallicity (right) at z = (top) and z = 2 (bottom). The colors and line styles are identical to those used in Fig. lCll 
The vertical, dotted lines in the left panels indicate the threshold for star formation. Note that star-forming gas was excluded from the temperature panels. The 
highest resolution simulations have nearly converged. 
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